The Magnus expansion and some of its 

applications 



S. Blanes% F. Casas^, J. A. Oteo'^ and J. Ros^'^ 

^Instituto de Matemdtica Multidisciplinar, Universidad Politecnica de Valencia, 

E-46022 Valencia, Spain 

^Departament de Matematiques, Universitat Jaume I, E-12071 Castellon, Spain 

Departament de Fisica Tedrica, Universitat de Valencia, E-46100 Burjassot, 

Valencia, Spain 

'^IFIC, Centre Mixt Universitat de Valencia- CSIC, E-46100 Burjassot, Valencia, 

Spain 



Abstract 

Approximate resolution of linear systems of differential equations with varying 
coefficients is a recurrent problem shared by a number of scientific and engineering 
areas, ranging from Quantum Mechanics to Control Theory. When formulated in 
operator or matrix form, the Magnus expansion furnishes an elegant setting to built 
up approximate exponential representations of the solution of the system. It provides 
a power series expansion for the corresponding exponent and is sometimes referred to 
as Time-Dependent Exponential Perturbation Theory. Every Magnus approximant 
corresponds in Perturbation Theory to a partial re-summation of infinite terms with 
the important additional property of preserving at any order certain symmetries of 
the exact solution. 

The goal of this review is threefold. First, to collect a number of developments 
scattered through half a century of scientific literature on Magnus expansion. They 
concern the methods for the generation of terms in the expansion, estimates of the 
radius of convergence of the series, generalizations and related non-perturbative 
expansions. Second, to provide a bridge with its implementation as generator of 
especial purpose numerical integration methods, a field of intense activity during 
the last decade. Third, to illustrate with examples the kind of results one can expect 
from Magnus expansion in comparison with those from both perturbative schemes 
and standard numerical integrators. We buttress this issue with a revision of the 
wide range of physical applications found by Magnus expansion in the literature. 
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1 Introduction 



1.1 Motivation, overview and history 



The outstanding mathematician Wilhelm Magnus (1907-1990) did important 
contributions to a wide variety of fields in mathematics and mathematical 
physics [1]. Among them one can mention combinatorial group theory [159] 
and his collaboration in the Bateman project on higher transcendental func- 
tions and integral transforms [81]. In this report we review another of his 
long-lasting constructions: the so-called Magnus expansion (hereafter referred 
to as ME). ME was introduced as a tool to solve non-autonomous linear dif- 
ferential equations for linear operators. It is interesting to observe that in his 
seminal paper of 1954 [158], although it is essentially mathematical in nature, 
Magnus recognizes that his work was stimulated by results of K.O. Friedrichs 
on the theory of linear operators in Quantum Mechanics [88]. Furthermore, 
as the first antecedent of his proposal he quotes a paper by R.P. Feynman 
in the Physical Review [86]. We stress these facts to show that already in its 
very beginning ME was strongly related to Physics and so has been ever since 
and there is no reason to doubt that it will continue to be. This is the first 
motivation to offer here a review as the present one. 

Magnus proposal has the very attractive property of leading to approximate 
solutions which exhibit at any order of approximation some qualitative or 
physical characteristics which first principles guarantee for the exact (but un- 
known) solution of the problem. Important physical examples are the symplec- 
tic or unitary character of the evolution operator when dealing with classical 
or quantum mechanical problems, respectively. This is at variance with most 
standard perturbation theories and is apparent when formulated in a correct 
algebraic setting: Lie algebras and Lie groups. But this great advantage has 
been some times darkened in the past by the difficulties both in constructing 
explicitly higher order terms and in assuring existence and convergence of the 
expansion. 

In our opinion, recent years have witnessed great improvement in this situa- 
tion. Concerning general questions of existence and convergence new results 
have appeared. From the point of view of applications some new approaches 
in old fields have been published while completely new and promising avenues 
have been open by the use of Magnus expansion in Numerical Analysis. It 
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seems reasonable to expect fruitful cross fertilization between these new de- 
velopments and the most conventional perturbative approach to ME and from 
it further applications and new calculations. 

This new scenario makes desirable for the Physics community in different areas 
(and scientists and engineers in general) to have access in as unified a way as 
possible to all the information concerning ME which so far has been treated 
in very different settings and has appeared scattered through very different 
bibliographic sources. 

As implied by the preceding paragraphs this report is mainly addressed to 
a Physics audience, or closed neighbors, and consequently we shall keep the 
treatment of its mathematical aspects within reasonable limits and refer the 
reader to more detailed literature where necessary. By the same token the 
applications presented will be limited to examples from Physics or from the 
closely related field of Physical Chemistry. We shall also emphasize its instru- 
mental character for numerically solving physical problems. 

In the present section as an introduction we present a brief overview and 
sketchy history of more than 50 years of ME. To start with, let us consider the 
initial value problem associated with the linear ordinary differential equation 

Y'{t)^A{t)Y{t), Y{to)^Yo, (1) 

where as usual the prime denotes derivative with respect to the real inde- 
pendent variable which we take as time t, although much of what will be 
said applies also for a complex independent variable. In order of increasing 
complexity we may consider the equation above in different contexts: 

(a) y : R — > C, ^ : R — > C. This means that the unknown Y and the 
given A are complex scalar valued functions of one real variable. In this 
case there is no problem at all: the solution reduces to a quadrature and 
an ordinary exponential evaluation: 

Y{t)=exp(^J^' A{s)ds^ Fo- (2) 

(b) Y :R — > C", A : M — ^ M„(C), where M„(C) is the set of n x n complex 
matrices. Now F is a complex vector valued function and A a complex 
n xn matrix valued function. At variance with the previous case, only in 
very special cases is the solution easy to state: when for any pair of values 
of t, ti and ^2, one has ^1(^1)74(^2) = ^(^2)^(^1), which is certainly the 
case if A is constant. Then the solution reduces to a quadrature (trivial or 
not) and a matrix exponential. With the obvious changes in the meaning 
of the symbols, equation (2) still applies. In the general case, however, 
there is not compact expression for the solution and (2) is not any more 
the solution. 
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(c) y : R — > Mn{C), A : R — ^ M„(C). Now both Y and A are com- 
plex matrix valued functions. A particular case, but still general enough 
to encompass the most interesting physical and mathematical applica- 
tions, corresponds to Y{t) G Q, A{t) G g, where Q and g are respectively 
a matrix Lie group and its corresponding Lie algebra. Why this is of 
interest is easy to grasp: the key reason for the failure of (2) is the non- 
commutativity of matrices in general. So one can expect that the (in 
general non-vanishing) commutators play an important role. But when 
commutators enter the play one immediately thinks in Lie structures. Fur- 
thermore, plainly speaking, the way from a Lie algebra to its Lie group 
is covered by the exponential operation. A fact that will be no surprise 
in this context. The same comments of the previous case are valid here. 
In this report we shall mostly deal with this matrix case. 

(d) The most general situation one can think of corresponds to Y{t) and A(t) 
being operators in some space, e.g., Hilbert space in Quantum Mechanics. 
Perhaps the most paradigmatic example of (1) in this setting is the time- 
dependent Schrodinger equation. 

Observe that case (b) above can be reduced to case (c). This is easily seen 
if one introduces what in mathematical literature is named the matrizant, a 
concept dating back at least to the beginning of the 20th century in the work 
of Baker [8]. It is the n x n matrix U{t,to) defined through 

Y{t)^U{t,to)Yo. (3) 

Without loss of generality we will take to = unless otherwise explicitly stated, 
for the sake of simplicity. When no confusion may arise, we write only one 
argument in U and denote U{t,0) = U{t), which then satisfy the differential 
equation and initial condition 

U'{t) = A{t)U{t), U{0)=I, (4) 

where / stands for the n-dimensional identity matrix. The reader will have 
recognized U (t) as what in physical terms is known as time evolution operator. 

We are now ready to state Magnus proposal: a solution to (4) which is a true 
matrix exponential 

u{t) = expn{t), n{o) = o, (5) 

and a series expansion for the matrix in the exponent 

oo 

n{t) = Y^nkit), (6) 

k=l 

which is what we call the Magnus expansion. The mathematical elaborations 
explained in the next section determine ilfe(t). Here we just write down the 
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three first terms of that series: 

Jo 

Mt) dti r di2 [A{tr),A{t2)] (7) 

/ JO Jo 

ns{t) = lfdt, f'dh t'dts {[A{U),[A{h),A{ts)]] + [A{ts),[A{t2),A{t,m 

D Jo Jo Jo 

where [A, B] = AB — BA is the matrix commutator of A and B. 

The interpretation of these equations seems clear: Qi{t) coincides exactly with 
the exponent in (2). But this equation cannot give the whole solution as has 
already been said. So, if one insists in having an exponential solution the 
exponent has to be corrected. The rest of the ME in (6) gives that correction 
necessary to keep the exponential form of the solution. 

The terms appearing in (7) already suggest the most appealing characteristic 
of ME. Remember that a matrix Lie algebra is a linear space in which one 
has defined the commutator as the second internal composition law. If, as we 
suppose, A{t) belongs to a Lie algebra q for all t so does any sum of multiple 
integrals of nested commutators. Then, if all terms in ME have a structure 
similar to that of the ones shown before, the whole fl{t) and any approximation 
to it obtained by truncation of ME will also belong to the same Lie algebra. 
In the next section it will be shown that this turns out to be the case. And a 
fortiori its exponential will be in the corresponding Lie group. 

Why is this so important for physical applications? Just because many of the 
properties of evolution operators derived from first principles are linked to 
the fact that they belong to a certain Lie group: e.g. unitary group in Quan- 
tum Mechanics, symplcctic group in Classical Mechanics. In that way use 
of (truncated) ME leads to approximations which share with the exact solu- 
tion of equation (4) important qualitative (very often, geometric) properties. 
For instance, in Quantum Mechanics every approximant preserves probability 
conservation. 

From the present point of view we could say that the last paragraphs summa- 
rize, in a nut shell, the main contents of the famous paper of Magnus of 1954. 
With no exaggeration its appearance can be considered a turning point in the 
treatment of the initial value problem defined by (4) . 

But important, as it certainly was, Magnus paper left some problems, at least 
partially, open: 

• First, for what values of t and for what operators A does equation (4) admit 
a true exponential solution? This we call the existence problem. 

• Second, for what values of t and for what hnear operators A does the series 
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in equation (6) converge? This we describe as the convergence problem. We 
want to emphasize that, although related, these two are different problems. 
To see why, think of the scalar equation y' = y"^ with y(0) = 1. Its solution 
y{t) = (1 — t)~^ exists for t ^ 1, but its form in power series y{t) = YIq" ^" 
converges only for \t\ < 1. 

• Third, how to construct higher order terms fife(t). A; > 3, in the series? 
Moreover, is there a closed- form expression for Vtk{t)l 

• Fourth, how to calculate in an efficient way exp VL^^\ where Vt^^^ = Yl,k=i ^k{t) 
is a truncation of the ME? 

All these questions, and many others, will be dealt with in the rest of the pa- 
per. But before entering that analysis we think interesting to present a view, 
however brief, from a historical perspective of this half a century of devel- 
opments on Magnus series. Needless to say, we by no means try to present a 
detailed and exhaustive chronological account of the many approaches followed 
by authors from very different disciplines. To minimize duplication with later 
sections we simply mention some representative samples in order the reader 
can understand the evolution of the field. 

Including some precedents, and with a, as undeniable as unavoidable, dose of 
arbitrariness, we may distinguish four periods in the history of our topic: 

(1) Before 1953. The problem which ME solves has a centennial history dat- 
ing back at least to the work of Peano, by the end of 19th century, and 
Baker, at the beginning of the 20th (for references to the original papers 
see e.g. [110]). They combine the theory of differential equations with an 
algebraic formulation. Intimately related to these treatments from the 
very beginning is the study of the so called Baker-Campbell-HausdorfT 
or BCH formula for short [9,40,102] which gives C in terms of A, B 
and their multiply nested commutators when expressing exp{A) exp(5) 
as exp(C). This topic has by itself a long history up until today and will 
be discussed also in section 2. As one of its early hallmarks we quote 
[75]. In Physics literature the interest in the problem posed by equation 
(4) highly revived with the advent of Quantum Electrodynamics (QED). 
The works of F.J. Dyson [76] and in particular R.P. Feynman [86] in the 
late forties and early fifties are worth mentioning here. 

(2) 1953-1970. We have quoted as birth certificate of ME the paper [158] by 
Magnus in 1954. This is not strictly true: there is a Research Report [157] 
dated June 1953 which differs from the published paper in the title and 
in a few minor details and which should in fact be taken as a preliminary 
draft of it. In both publications appears the result summarized above 
on ME with almost identical words. The work of Pechukas and Light 
[195] gave for the first time a more specific analysis of the problem of 
convergence than the rather vague considerations in Magnus paper. Wei 
and Norman [235,236] did the same for the existence problem. Robinson, 
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to the best of our knowledge, seems to have been the first to apply ME 
to a physical problem [203]. Special mention in this period deserves a 
paper by Wilcox [239] , in which useful mathematical tools are given and 
ME presented together with other algebraic treatments of equation (4), 
in particular Fer's infinite product expansion [85] . Worth also of mention 
here is the first application of ME as a numerical tool for integrating the 
time-independent Schrodinger equation for potential scattering by Chang 
and Light [50]. 

(3) 1971-1990. During these years ME consolidated in different fronts. It was 
successfully apphed to a wide spectrum of fields in Physics and Chem- 
istry: from atomic [11] and molecular [210] Physics to Nuclear Magnetic 
Resonance (NMR) [82,234] to Quantum Electrodynamics [59] and ele- 
mentary particle Physics [68]. A number of case studies also helped to 
clarify its mathematical structure, see for example [135]. The construc- 
tion of higher order terms was approached from different angles. The 
intrinsic and growing complexity of ME allows for different schemes. One 
which has shown itself very useful in tackling other questions like the 
convergence problem was the recurrent scheme by Klarsfeld and Oteo 
[137]. 

(4) Since 1991. The last decade of the 20th century witnessed a renewed inter- 
est in ME which still continues nowadays. It has followed different lines. 
Concerning the basic problems of existence and convergence of ME there 
has been definite progress [21,43,173,176]. ME has also been adapted for 
specific types of equations: Floquet theory when A{t) is a periodic func- 
tion [45], stochastic differential equations [34] or equations of the form 
Z' = AZ — ZB [112]. Special mention should be made to the new field 
open in this most recent period that uses Magnus scheme to build novel 
algorithms [120] for the numerical integration of differential equations 
within the most wide field of geometric integration [32]. After optimiza- 
tion [22,23], these integrators have proved to be highly competitive. 

As a proof of the persistent impact the 1954 paper by Magnus has had in 
scientific literature we present in Figures 1 and 2 the number of citations per 
year and the cumulative number of citations, respectively, as December 2007 
with data taken from ISI Web of Science. The original paper appears about 
750 times of which, roughly, 50, 320 and 380 correspond respectively to each 
of the last three periods we have considered. The enduring interest in that 
seminal paper is clear from the figures. 

The presentation of this report is organized as follows. In the remaining of 
this section we include some mathematical tools and notations that will be 

used time and again in our treatment. In section 2 we introduce formally 
the Magnus expansion, study its main features and analyze thoroughly the 
convergence issue. Next, in section 3 several generalizations of the Magnus 
expansion are reviewed, with special emphasis in its application to general 
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Fig. 1. Persistency of Magnus original paper: number of citations per year. 
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Fig. 2. Persistency of Magnus original paper: cumulative number of citations. 

nonlinear differential equations. In order to illustrate the main properties of 
ME, in section 4 we consider simple examples for which the computations 
required are relatively straightforward. Section 5 is devoted to an aspect that 
has been most recently studied in this setting: the design of new algorithms 
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for the numerical integration of differential equations based on the Magnus 
expansion. There, after a brief characterization of numerical integrators, we 
present several methods that are particularly efficient, as shown by the exam- 
ples considered. Given the relevance of the new numerical schemes, we briefly 
review in section 6 some of its applications in different contexts, ranging from 
boundary-value problems to stochastic differential equations. In section 7, on 
the other hand, applications of the ME to significant physical problems are 
considered. Finally, the paper ends with some concluding remarks. 

1.2 Mathematical preliminaries and notations 

Here we collect for the reader's convenience some mathematical expressions, 
terminology and notations which appear most frequently in the text. Needless 
to say that we have made no attempt to being completely rigorous. We just 
try to facilitate the casual reading of isolated sections. 

As already mentioned, the natural mathematical habitat for most of the ob- 
jects we will deal with in this report is a Lie group or its associated Lie 
algebra. Although most of the results discussed in these pages are valid in a 
more general setting we will essentially consider only matrix Lie groups and 
algebras. 

By a Lie group Q we understand a set which combines an algebraic structure 
with a topological one. At the algebraic level every two elements of Q can be 
combined by an internal composition law to produce a third element also in 
Q. The law is required to be associative, to have an identity element and every 
element must have an inverse. The ordinary product and inverse of invertible 
matrix play that role in the cases we are more interested in. The topological 
exigence forces the composition law and the association of an inverse to be 
sufficiently smooth functions. 

A Lie algebra g is a vector space whose elements can be combined by a second 
law, the Lie bracket, which we represent by [A, B] = C, with A, B, C elements 
of 0, in such a way that the law is bilinear, skew-symmetric and satisfies the 
well known Jacobi identity, 

[A, [B, C]] + [B, [C, A]] + [C, [A, B]] = 0. (8) 

When dealing with matrices we take as Lie bracket the familiar commutator: 

[A,B\^AB-BA, AeQ, B e Q, (9) 

where AB stands for the usual matrix product. If we consider a finite-dimen- 
sional Lie algebra with dimension d and denote hy Ai,i = 1, ... ,d, the vectors 
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of one of its basis then the fundamental brackets one has to know are 



[A, A,] = cf^.Afc, (10) 

where sum over repeated indexes is understood. The coefficients c^j are the 
so-called structure constants of the algebra. 

Associated with any A & Qwe can define a linear operator ad^ '■ 9 ^ Q which 
acts according to 

ad^S = [A, B], ad^S = [A, ad^"^S], ad^B ^B, j eN,B eg. 

(11) 

Also of interest is the exponential of this ad^ operator, 

AdA = exp(adA), (12) 

whose action on q is given by 

OO 

AdA{B) = exp(A) B exp(-A) = ^ -ryad^S, Beg. (13) 

fe=0 ^' 

The type of matrices we will handle more frequently are orthogonal, unitary 
and symplectic. Here are their characterization and the notation we shall use 
for their group and algebra. 

The special orthogonal group, SO(n), is the set of all n x n real matrices with 
unit determinant satisfying A^A = AA^ = I, where A^ is the transpose of A 
and / denotes the identity matrix. The corresponding algebra so{n) consists 
of the skew-symmetric matrices. 

A n X n complex matrix A is called unitary if A^A = AA^ = I, where A^ is 
the conjugate transpose or Hermitian adjoint of A. The special unitary group, 
SU(n), is the set of all n x n unitary matrices with unit determinant. The 
corresponding algebra su(n) consist of the skew-Hermitian traceless matrices. 
Special relevance in some quantum mechanical problems we discuss will have 
the case n = 2. In this case a convenient basis for 5u(2) is made up by the 
Pauli matrices 



^2 = , ^3 = . (14) 






They satisfy the identity 

ajak = Sjk + iejkiCJi, (15) 

and correspondingly 

[aj,ak] = 2iejkiai, (16) 
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which directly give the structure constants for SU(2). The following identities 
will prove useful for a and 6 in M^: 

{a ■ cr){b ■ cr) ^ a - b I + i{a x b) ■ cr, [a • cr, 6 • <t] = 2i{a x b) ■ (t, (17) 

where we have denoted cr = (cxi, (T2, a^). Any U G SU(2) can be written as 

/• \ / \ T .sin(a) 
U — expna ■ cr) — cos(a) / + z — —a ■ cr, (18) 

a 

where a = \\a\\ ~ ^Ja^~+~a^'+~a^. A more elaborate expression which we shall 
make use of in later sections is (with a — 1) 

exp(ia-crt) (b-a) exp(—ia-crt) — b-cr+sm2t {bxa)-cr+sm'^ t {ax{bxa))-cr 

(19) 

In Hamiltonian problems the symplectic group Sp(n) plays a fundamental role. 
It is the group of 2n x 2n real matrices satisfying 



A^JA = J, with J = 



^ On In 



'In 



(20) 



and /„ denotes the n-dimensional identity matrix. Its corresponding Lie alge- 
bra sp(n) consists of matrices verifying B^J + JB = 02n- In fact, these can be 
considered particular instances of the so-called J-orthogonal group, defined as 
[197] 

Oj{n) = {A e GL(n) : A^JA = J}, (21) 
where GL{n) is the group of all n x n nonsingular real matrices and J is some 
constant matrix in GL(n). Thus, one recovers the orthogonal group when 
J = I, the symplectic group Sp(n) when J is the basic symplectic matrix 
given in (20), and the Lorentz group S0(3, 1) when J = diag(l, —1, —1, — 1). 
The corresponding Lie algebra is the set 

oj{n) = {Be 0[„(M) : B^J + JB = O}, (22) 

where g[„(M) is the Lie algebra of all n x n real matrices. If S e Oj(n), then 
its Cayley transform 

A = {I - aB)-\l + aB) (23) 

is J-orthogonal. 

Another important matrix Lie group not included in the previous characteri- 
zation is the special linear group SL(n), formed by all n x n real matrices with 
unit determinant. The corresponding Lie algebra sl(n) comprises all traceless 
matrices. For real 2x2 matrices in s 1(2) one has 



exp 



( a 




\ c —a . 





cosh(77) -I- ^ sinh(77) ^ sinh(77) ^ 
^ sinh(77) cosh(77) - ^ sinh(77) ^ 



(24) 
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with T) — y/a^ + be. 



When deahng with convergence problems it is necessary to use some type 
of norm for a matrix. By such we mean a non- negative real number \\A\\ 
associated with each matrix A e C"^" and satisfying 

a) \\A\\ > for all A and P|| = iff A = 0„. 

b) \\q:A\\ = \a\ \\A\\, for all scalars a. 

c) \\A + B\\ < 11^11 + 

Quite often one adds the sub-multiplicative property 

\\AB\\<\\A\\\\Bl (25) 
but not all matrix norms satisfy this condition [93]. 

There exist different families of matrix norms. Among the more popular ones 
we have the p-norm \\A\\p and the Frobenius norm For a matrix A with 

elements a^j, i,j — 1 . . .n, they are defined as 



ll^llp^max II^xIIp (26) 

X L=l 



n n 



i=i j=i 



respectively, where ||x||p = {J2]=i kj|^)^ and tr(A) is the trace of the matrix 
A. Although both verify (25), the p-norms have the important property that 
for every matrix A and x G R" one has ||^x||p < \\A\\p ||x||p. The most used 
p-norms correspond to p = 1, p = 2 and p — oo. 

Of paramount importance in numerical linear algebra is the case p — 2. The 
resulting 2-norm of a vector is nothing but the Euclidean norm, whereas in the 
matrix case it is also called the spectral norm of A and can be characterized as 
the square root of the largest eigenvalue of A^A. A frequently used inequality 
relating Frobenius and spectral norms is the following: 

\\A\U<\\A\W<V^\\A\U. (28) 

In fact, this last inequality can be made more stringent [227]: 



\\A\\f < Jrank(A) \\A\\2. (29) 



Considering in a matrix Lie algebra q a norm satisfying property (25), it is clear 
that ||[74,5]|| < 2||yl||||i?||, and the ad operator defined by (11) is bounded. 



14 



since 

llad^ll < 2||A|| 

for any matrix A. 

A matrix norm is said to be unitarily invariant if ||J7Ay|| = \\A\\ whenever U, 
V are unitary matrices. Probenius and p-norms are both unitarily invariant 
[107]. 

In some of the most basic formulas for the Magnus expansion there will appear 
the so-called Bernoulli numbers B^, which are defined through the generating 
function [2] 

^^EBniz)'-., \t\<2n 

c* — 1 nl 



n=0 

as Bn = Bn{0). Equivalent ly, 



r °° B 

e--l to ^! 

whereas the formula 

f^x 1 oo 1 

! i ^ V = 

^ ^o(^ + l)! 

will be also useful in the sequel. The first few nonzero Bernoulli numbers are 
Bq = 1, Bi = — i, B2 = |, 54 = — ^. In general one has -82^+1 = for m > 1. 



2 The Magnus expansion (ME) 



Magnus proposal with respect to the linear evolution equation 

Y'{t) = A{t)Y{t) (30) 

with initial condition 1^(0) = /, was to express the solution as the exponential 
of a certain function, 

y(i) = expQ(i). (31) 
This is in contrast to the representation 

Y{t) = T (exp A{s)ds^ 

in terms of the time-ordening operator T introduced by Dyson [76]. 

It turns out that Q{t) in (31) can be obtained explicitly in a number of ways. 
The crucial point is to derive a differential equation for the operator Q that 
replaces (30). Here we reproduce the result first established by Magnus as 
Theorem III in [158]: 
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Theorem 1 (Magnus 1954)- Let A{t) he a known function oft (in general, in 
an associative ring), and let Y{t) he an unknown function satisfying (30) with 
F(0) = /. Then, if certain unspecified conditions of convergence are satisfied, 
Y{t) can he written in the form 

Y{t) = expO(i), 

where 

^ = V^ad^A, (32) 
dt t,n\ 

and Bn are the Bernoulli numbers. Integration of (32) hy iteration leads to an 
infinite series for fl the first terms of which are 



n{t) = f A{ti)dti - - / 

Jo 2 Jo 

2.1 A proof of Magnus Theorem 



A{t2)dt2,A{ti) 







dti + 



The proof of this theorem is largely based on the derivative of the matrix 
exponential map, which we discuss next. Given a scalar function uj{t) G M, 
the derivative of the exponential is given by dexp{uj(t))/dt = uj'it) exp(a;(t)). 
One could think of a similar formula for a matrix VL{t). However, this is not 
the case, since in general [Q, Q!] ^ 0. Instead one has the following result. 

Lemma 2 The derivative of a matrix exponential can he written alternatively 
as 



(a) ^exp{Q{t)) = dexp^^,^{Q'{t))exp{Q{t)), (33) 



dt 

(b) j^exp{Q{t)) = exp{Q{t))dexp_^^,){Q'{t)), (34) 

(c) 4- exp(0(i)) = r e^^W n'{t) e(i-^)"W dx, (35) 
dt Jo 

where dexp^{C) is defined hy its (everywhere convergent) power series 

dexMC) = E 7r^ad^(C) ^ ^M^y-I ^c). (36) 

Proof. Let 0,{t) be a matrix- valued differentiable function and set 

d 

Y{a,t) = — {exp{an{t)))exp{-an{t)) 
for (7, i e R. Differentiating with respect to cr. 
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dY d d 

= — {cxp{aQ)Q) exp(— (7^2) + — (cxp(crr2)) (— fi) cxp(— 

f d \ 

— I exp((7Q)Q' + — (exp((TQ)) Q j exp(— (tQ) 

d 

— — (exp((7n)) Oexp(— crO) = exp{ail)il' exp{—ail) 

oo fc 

= exp(ad.o)(l^') = E7T^4(^'), 
fe=o ■ 

where the first equahty in the last line follows readily from (12) and (13). On 
the other hand 

4(expQ)exp(-Q) = y(l,i) = C ■^Y(a,t)da (37) 
dt Jo oa 

since Y{0, t) — 0, and 



from which formula (33) follows. The convergence of the power series (36) is 
a consequence of the boundedness of the ad operator: Hadnll < 2||0||. 

Multiplying both sides of (33) by exp(— O), we have 

e-" —- = e-" dexpj^(n') e" = e^^^-" d&cp^{n') = Q' = d&cp_^{n') 

dt ad_n 



from which (34) follows readily. Finally, equation (35) is obtained by taking 

-Y(a,t)da — [ exp(aD,)D,' ex.p(—aD,)da 
Jo 



da 



in (37). 



According to Rossmann [205] and Sternberg [218], formula (33) was first 
proved by F. Schur in 1890 [212] and was taken up later from a different 
point of view by Poincare (1899), whereas the integral formulation (35) has 
been derived a number of times in the physics literature [239] . 

As a consequence of the Inverse Function Theorem, the exponential map has a 
local inverse in the vicinity of a point at which d exp^ = (exp(adf2 ) — /)/ adf^ 
is invertible. The following lemma establishes when this takes place. 

Lemma 3 (Baker 1905). If the eigenvalues of the linear operator ad^ are 
different from 2mm with m e {±1,±2,...}, then dexp^ is invertible. Fur- 
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thermore, 

^-P^^(^) = = E ^-^-(C) (38) 

and the convergence of the dexp^^ expansion is certainly assured if \\^\\ < tt. 
Proof. The eigenvalues of d exp^ are of the form 

where v is an eigenvalue of adf^. By assumption, the values of ^ are non- 
zero, so that (iexpfj is invertible. By definition of the Bernoulli numbers, the 
composition of (38) with (36) gives the identity. Convergence for \\VL\\ < tt 
follows from Hadf^H < 2||Q|| and from the fact that the radius of convergence 
of the series expansion for x/(e^ — 1) is 27r. ■ 

It remains to determine the eigenvalues of the operator ad^- In fact, it is not 
difficult to show that if Q has n eigenvalues {Xj, j — 1,2, ... ,n}, then adf^ 
has eigenvalues {Xj — Xk, j,k — 1,2, ... ,n}. 

As a consequence of the previous discussion. Theorem 1 can be rephrased more 
precisely in the following terms. 

Theorem 4 The solution of the differential equation Y' = A(t)Y with initial 
condition Y{0) = Yq can be written as Y{t) — exp{fl{t))Yo with fl{t) defined 
by 

Q' = dexp^\A{t)), Q{0) = O, (39) 

where 

oo r) 

rfexp^^(>l) = ^-fad^(A). 

k=Q ■ 

Proof. Comparing the derivative of Y{t) — exp{Q,{t))Yo, 

^ = ^ (exp(Q(t)))yo = ^^exp^(Q') eMm)yo 

with Y' = A(t)Y, we obtain A(t) = d expQ[Q') . Applying the inverse operator 
dexp^^ to this relation yields the differential equation (39) for fl{t). m 

Taking into account the numerical values of the first few Bernoulli numbers, 
the differential equation (39) therefore becomes 

n' = A(t) - ^[n, A(t)] + ^[n, [n, A{t)]] + ■■■, 
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which is nonhnear in Q. By defining 

= O, = f' A(tMti, 

Jo 

and applying Picard fixed point iteration, one gets 

and hm„_^oo ^'"'(^) = ^{t) in a suitably small neighbourhood of the origin. 
2.2 Formulae for the first terms in Magnus expansion 

Suppose now that A is of first order in some parameter e and try a solution 
in the form of a series 

oo 

0(t) = ^0„(t), (40) 

n=l 

where 0„ is supposed to be of order e^. Equivalently, we replace A i — > sA in 
(30) and determine the successive terms of 

oo 

0(t) = ^£"0„(t). (41) 

n=l 

This can be done cxphcitly, at least for the first terms, by substituting the 
series (41) in (39) and equating powers of £. Obviously, the Magnus series (40) 
is recovered by taking e — 1. Thus, using the notation A{ti) = Ai, the first 
four orders read 

(1) = A, so that 

VLAt) = I dtiAi (42) 

(2) = Thus 

n2{t)^l fdti rdt2[A,,A2] (43) 
2 Jo Jo 

(3) = —^[fl2, A] + -^[fli, [fli, A]]. After some work and using the formula 

rOi fX ra ra 

dx f{x,y)dy= dy f{x,y)dx (44) 

Jo Jo Jo Jy 

we obtain 

n^{t) = \ fdt, tdt^ rdts{[A,,[A2,As]] + [[A,,A2lAs]} (45) 

D Jo Jo Jo 
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(4) = -|[Q3, A\ + ^[^2, [111, A]] + ii[^li, [^^2, A\\, which yields 

04(i) ^^f^ dii dt2 dh di4{[[[^l, ^2], A3]>l4] (46) 
+ [Ai, A3], A4]] + [Ai, [X, [A3, A4]]] + [A2, [A3, [A4, Ai]]]} 

The apparent symmetry in the formulae above is deceptive. High orders require 
repeated use of (44) and become unwieldy. Prato and Lamberti [198] give 
explicitly the fifth order using an algorithmic point of view. One can also 
find in the literature quite involved explicit expressions for an arbitrary order 
[16,167,206,219,220]. In the next subsection we describe a recursive procedure 
to generate the terms in the expansion. 



2.3 Magnus expansion generator 



The above procedure can provide indeed a recursive procedure to generate all 

the terms in the Magnus series (40). Thus, by substituting Vt{t) = 

into equation (39) and equating terms of the same order one gets in general 

fl[ = A 

K = j:^S^n'\ n>2, (47) 
j=i J- 

where 

Si''> = Y.[n^„[■■■[^^,,A]...]] (^, + ... + ^, = r^-l). (48) 

Notice that in the last equation the order in A has been explicitly reckoned, 
whereas k represents the number of Jl's. The newly defined operators Sj^^ can 
again be calculated recursively. The recurrence relations are now given by 



Si'^ ^ E [^m, Strlr^] , 2<j<n-l (49) 

m=l 

5« = [Q„_i,A], 5M^ad^7^(A). 
After integration we reach the final result in the form 



Qi= /*A(r)d7 
Jo 



'^^^1 R ft 

^n=E-f / S^n\r)dT, n>2. (50) 



j\ Jo 
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Alternatively, the expression of S^^ given by (48) can be inserted into (50), 
thus arriving at 

^n{t) = Y.-f Yl / adn,^(.) adn,^(,) • • • adn,.(,)^(s) ds n>2. 

j=l J- fcjH \-kj=n-l "'O 

ki>l,...,kj>l 

(51) 

Notice that each term n„(i) in the Magnus series is a multiple integral of 
combinations of n — 1 nested commutators containing n operators A{t). If, 
in particular, A{t) belongs to some Lie algebra 0, then it is clear that Q{t) 
(and in fact any truncation of the Magnus series) also stays in q and therefore 
exp(r2) G Q, where Q denotes the Lie group whose corresponding Lie algebra 
(the tangent space at the identity of Q) is q. 



2.4 Magnus expansion and time- dependent perturbation theory 

It is not difficult to establish a connection between Magnus series and Dyson 
perturbative series [76]. The later gives the solution of (30) as 

00 

Y{t)^I+^Pn{t), (52) 

n=l 

where P„ are time-ordered products 

Pn{t) = f^dti... ' dtn A1A2 ...An, 

Jo Jo 

where A = A{ti). Then 

00 / 00 

5:^(i)=iog /+^p,(t) 

i=i V j=i , 

As stated by Salzman [208], 

nn = Pn-Y.—^Rn\ ^ > 2, (53) 

j=2 J 

where 

obeys to the quadratic recursion formula 



rI^^^^'Er^^^rH:^, (54) 

m=l 

d(1) _ p r>[n) _ pn 
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Equation (54) represents the Magnus expansion generator in Salzman's ap- 
proach. It may be useful to write down the first few equations provided by 
this formahsm: 



^2^P2-\PI (55) 

A similar set of equations was developed by Burum [36], thus providing 
Pi^^i, 

P2 = Q2 + ^^^?, (56) 

P3 = ^3 + ^(OiOs + ^2^1) + J^^l 

and so on. The general term reads 

nn^Pn-Y.-Q^\ n>% (57) 

where 

= 5^ Qi, . . . (ii + • • • + ife = n). (58) 

As before, subscripts indicate the order with respect to the parameter £, 
while superscripts represent the number of factors in each product. Thus, 
the summation in (58) extends over all possible products of k (in general non- 
commuting) operators Jlj such that the overall order of each term is equal to 
n. By regrouping terms, one has 

i2-\ Hk=n-l i2-\ \-ik=n-'2 

-|- • • • -|- ^ ^ ■ ■ ■ 1 

i2-\ |-jfe=fc— 1 

where Q^-* niay also be obtained recursively from 



n-j+l 

Q^^= E Q^^Qt'nl, (60) 

m=l 

Qn ^ — Qn ^ — 



By working out this recurrence one gets the same expressions as (55) for the 
first terms. Further aspects of the relationship between Magnus, Dyson series 
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and time-ordered products can be found in [143] and [192]. 



2.5 Graph theoretical analysis of Magnus expansion 



The previous recursions allow us in principle to express any in the Magnus 
series in terms of Qi, . . . ,Qk~i- In fact, this procedure has some advantages 
from a computational point of view. On the other hand, as we have mentioned 
before, when the recursions are solved explicitly, flk can be expanded as a lin- 
ear combination of terms that are composed from integrals and commutators 
acting iteratively on A. The actual expression, however, becomes increasingly 
complex with k, as it should be evident from the first terms (42)- (46). An 
alternative form of the Magnus expansion, amenable also for recursive deriva- 
tion by using graphical tools, can be obtained by associating each term in the 
expansion with a binary rooted tree, an approach worked out by Iserles and 
N0rsett [120]. For completeness, in the sequel we show the equivalence of the 
recurrence (49)- (50) with this graph theoretical approach. 

In essence, the idea of Iserles and N0rsett is to associate each term in 0,^ with 
a rooted tree, according to the following prescription. 

Let Tq be the set consisting of the single rooted tree with one vertex, then 
Tq = {•}, establish the relationship between this tree and A through the map 

• ^ A{t) 

and define recursively 



T 



Next, given two expansion terms Hj-^ and Hj^, which have been associated 
previously with ti G 7^.^ and T2 G T^j, respectively {ki + k2 — m — 1), we 
associate 



Hr{t) 



f H,^{i)di,H,,{t) 
Jo 



with 



7-2 



Thus, each H^- for r e 7^ involves exactly m integrals and m commutators. 



These composition rules establish a one-to-one relationship between a rooted 
tree t E T = Urn>o%n, and a matrix function i?T-(^) involving A, multivariate 
integrals and commutators. 
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Prom here it is easy to deduce that every r E %n, m > 1, can be written in a 
unique way as 

^3 I 

T2 




or T = a{Ti,T2, ... ,Ts). Then the Magnus expansion can be expressed in the 
form [119,120] 

OO f 

m = E E «(^) / HriOd^, (61) 

m=0 TGTm 

with the scalar q;( • ) = 1 and, in general, 

B * 

1=1 

Let us illustrate this procedure by writing down explicitly the first terms in 
the expansion in a tree formalism. In Ti we only have ki — k2 — 0, so that a 
single tree is possible. 



with q;(t) = —1/2. In 7^ there are two possibilities, namely ki = 0, k2 — 1 
and ki — 1, k2 — 0, and thus one gets 



Ti = • , T2 = \/ ^ T = \/ , q;(t) = ^ 
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Ti = T2 = • =^ T = \/ Q;(t) = I 



and the process can be repeated for any 7^. The correspondence between trees 
and expansion terms should be clear from the previous graphs. For instance, 
the last tree is nothing but the integral of A, commuted with A, integrated and 
commuted with A. In that way, by truncating the expansion (61) at m — 2 
we have 




+ •••, (62) 
i.e., the exphcit expressions collected in subsection 2.2. 
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Finally, the relationship between the tree formalism and the recurrence (49)- 
(50) can be established as follows. From (61) we can write 

oo R 

^ ^ a(r)i/.(t) = ^-f ^ a{T,)---a{Ts)H^^^,,...,r.y 

m=lTeTm s=l ■ ki,...,ks neTk. 

fciH \-ks=m—s 

Thus, by comparing (50) and (61) we have 

n^{t) = ^ a{r) / HriOdC = E -f / S^^^iOdC 

reTm-i ''^ 3=1 J- 

SO that 

•^m = E E «('^l)---«('0-)^a(n,...,r^)- 

ki,...,kj n&T^.. 
feiH \-kj=m-l-j 

In other words, each term 5"^^^ in the recurrence (49) carries on a complete 
set of binary trees. Although both procedures are equivalent, the use of (49) 
and (50) can be particularly well suited when high orders of the expansion are 
considered, for two reasons: (i) the enormous number of trees involved for large 
values of m and (ii) in (61) many terms are redundant, and a careful graph 
theoretical analysis is needed to deduce which terms have to be discarded 
[120]. 

Recently, an ME-type formahsm has been developed in the more abstract 
setting of dendriform algebras. This generalized expansion incorporates the 
usual one as a limit, but is formulated more in line with (non-commutative) 
Butcher series. In this context, the use of planar rooted trees to represent the 
expansion and the so-called pre-Lie product allows one to reduce the number 
of terms at each order in comparison with expression (61) [77]. 

2. 6 Time-symmetry of the expansion 

The map ip^ : F(to) — ^ ^(^) corresponding to the linear differential equation 
(30) with l^(to) = ^0 is time symmetric, (/?~* o y?* = Id, since integrating (30) 
from t = to to t = t f ioi every tf > to and back to to leads us to the original 
initial value Y{to) = Yq. Observe that, according with (3), the map can be 
expressed in terms of the fundamental matrix (or evolution operator) U {t, to) 
as ip*f{Yo) — U{tf,to)Yo. Then time-symmetry establishes that 

U{to,tf) = U-\tf,to) 
or, in terms of the Magnus expansion, 

n{tf,to) = -n{to,tf). 
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To take advantage of this feature, let us write the solution of (30) at the final 
time tf — to + s as 

yiti/2 + |) = exp (n{ty2 + \Mi2 - ^)) 1^(^1/2 - (63) 
where ti/2 = (to + t/)/2- Then 

5^(^1/2 - |) = exp (-Q(ii/2 + |, ti/2 - |)) Y{tii2 + |). (64) 
On the other hand, the solution at to can be written as 

Y{t^l2 - |) = exp (l](tv2 - \Mi2 + |)) Y{t^i2 + |), (65) 
so that, by comparing (64) and (65), 

^{txl2 - \MI2 + |) = -^{txl2 + |, tl/2 - |) (66) 

and thus ^2 docs not contain even powers of s. If Aif) is an analytic function 
and a Taylor series centered around ti/2 is considered, then each term in ^2^. 
is an odd function of s and, in particular, f22i+i(s) = 0(s^*+^). This fact has 
been noticed in [122,182] and will be fully exploited in section 5.4 when ana- 
lyzing the Magnus expansion as a numerical device for integrating differential 
equations. 

2. 7 Convergence of the Magnus expansion 

As we pointed out in the introduction, from a mathematical point of view, 
there are at least two different issues of paramount importance at the very 
basis of the Magnus expansion: 

(1) [Existence) For what values of t and for what operators A docs equation 
(30) admit an exponential solution in the form Y{t) = exp{Q{t)) for a 
certain Q{t)7 

(2) (Convergence) Given a certain operator A{t), for what values of t does 
the Magnus series (40) converge? In other words, when fl{t) in (31) can 
be obtained as the sum of the series (40)? 

Of course, given the relevance of the expansion, both problems have been 
extensively treated in the literature since Magnus proposed this formalism 
in 1954. We next review some of the most relevant contributions available re- 
garding both aspects, with special emphasis on the convergence of the Magnus 
series. 
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2.1.1 On the existence ofi}{t) 



In most cases one is interested in the case where A belongs to a Lie algebra 
g under the commutator product. In this general setting the Magnus theorem 
can be formulated as four statements concerning the solution of Y' — A{t)Y, 
each one more stringent than the preceding [235]. Specifically, 

(A) The differential equation Y' — A{t)Y has a solution of the form Y{t) ~ 
expQ(i). 

(B) The exponent Vt(t) lies in the Lie algebra q. 

(C) The exponent Vt{t) is a continuous differentiable function of A{t) and t, 
satisfying the nonlinear differential equation Q! = d exp^^ {A{t)) . 

(D) The operator n{t) can be computed by the Magnus series (40). 

Let us analyze in detail now the conditions under which statements (A)-(D) 



(A) If A{t) and Y{t) are n x n matrices, from well-known general theorems 
on differential equations it is clear that the initial value problem defined by 
(30) and Y(0) — I always has a uniquely determined solution Y{t) which is 
continuous and has a continuous first derivative in any interval in which A{t) 
is continuous [53] . Furthermore, the determinant of Y is always different from 
zero, since 



On the other hand, it is well known that any matrix Y can be written in the 
form expf2 if and only if detF ^ [91, p. 239], so that it is always possible 
to write Y{t) = expfl{t). 

In the general context of Lie groups and Lie algebras, it is indeed the regu- 
larity of the exponential map from the Lie algebra g to the Lie group Q that 
determines the global existence of an Q{t) G Q [66,207]: the exponential map 
of a complex Lie algebra is globally one-to-one if and only if the algebra is 
nilpotent, i.e. there exists a finite n such that adxiSidx2 • • • ^xn-i^n — 0, where 
X j are arbitrary elements from the Lie algebra. In general, however, the injec- 
tivity of the exponential map is only assured for ^ G such that ||^|| < pg for 
a real number pg > and some norm in q [173,174]. 

(B) Although in principle pg constitutes a sharp upper bound for the mere 
existence of the operator Q e g, its practical value in the case of differential 
equations is less clear. As we have noticed, any nonsingular matrix has a 
logarithm, but this logarithm might be in Ql{n, C) even when the matrix is 
real. The logarithm of Y{t) may be complex even for real A(t) [235]. In such a 
situation, the solution of (30) cannot be written as the exponential of a matrix 
belonging to the Lie algebra over the field of real numbers. One might argue 



hold. 




27 



that this is indeed possible over the field of complex numbers, but (i) the 
element Q cannot be computed by the Magnus series (40), since it contains 
only real rational coefficients, and (ii) examples exist where the logarithm of 
a complex matrix does not lie in the corresponding Lie subalgebra [235]. 

It is therefore interesting to determine for which range of t a real A{t) in (30) 
leads to a real logarithm. This issue has been tackled by Moan in [174] in the 
context of a complete normed (Banach) algebra, proving that if 



then the solution of (30) can be written indeed as Y{t) — expQ{t), where Q{t) 
is in the Banach algebra. 

(C) In his original paper [158], Magnus was well aware that if the function 
fl{t) is assumed to be differentiable, it may not exist everywhere. In fact, he 
related the differentiability issue to the problem of solving dexp^^fl') = A{t) 
with respect to il' and provided an implicit condition for an arbitrary A. 
More specifically, he proved the following result for the case oi n xn matrices 
(Theorem V in [158]). 

Theorem 5 The equation A{t) — dexpQ(Q') can be solved by Q' — dexp^^ A{t) 
for an arbitrary n x n matrix A if and only if none of the differences between 
any two of the eigenvalues ofD, equals 2Trim, where m — ±1, ±2, . . ., (m ^ 0). 

This result can be considered, in fact, as a reformulation of Lemma 3, but, 
unfortunately, has not very much practical application unless the eigenvalues 
of Q can easily be determined from those of A(t). One would like instead to 
have conditions based directly on A. 

2.7.2 Convergence of the Magnus series 

For dealing with the validity of statement (D) one has to analyze the conver- 
gence of the series J2T=i ^k- Magnus also considered the question of when 
the series terminates at some finite index m, thus giving a globally valid 
= Qi + • • • + Q^. This will happen, for instance, if 



identically for all values of t, since then fl^ = for k > 1. A sufficient (but 
not necessary) condition for the vanishing of all terms flk with k > n is that 




(67) 




[A{s,), [A{s2), [A{ss), ■ ■ ■ , [Ais^), A{sr.+,)] ■■•]]] = 
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for any choice of Si, . . . , Sn+i- In fact, the termination of the series cannot be 
estabhshed solely by consideration of the commutativity of A{t) with itself, 
and Magnus considered an example illustrating this point. 

In general, however, the Magnus series does not converge unless A is small in a 
suitable sense. Several bounds to the actual radius of convergence in terms of 
A have been obtained in the literature. Most of these results can be stated as 
follows. If ^m{t) denotes the homogeneous element with m — 1 commutators 
in the Magnus series as given by (51), then 0,{t) — J2m=i ^m{t) is absolutely 
convergent for < i < T, with 

r = max|i >0 : ^* ||A(s)||2 cis < r^j . (68) 

Thus, both Pechukas and Light [195] and Karasev and Mosolova [130] ob- 
tained Tc = log 2 — 0.693147..., whereas Chacon and Fomenko [48] got 
= 0.57745 .... In 1998, Blanes et al. [21] and Moan [172] obtained inde- 
pendently the improved bound 



= 1/ 

2 Jo 



2n I 

^^^^^-^dx. 5 = 1.08686870... (69) 



by analyzing the recurrence (49)- (50) and (51), respectively. Furthermore, 
Moan also obtained a bound on the individual terms of the Magnus series 
[174] which is useful, in particular, for estimating errors when the series is 
truncated. Specifically, he showed that 

where fm are the coefficients of 

G~^(x) = y fmX"^ = x+-X^ + —X^ + —X'^ + -^^X^ -I , 

4 72 576 86400 

the inverse function of 

rs 1 

dx. 



cot I) 



On the other hand, by analyzing some selected examples. Moan [174] con- 
cluded that, in order to get convergence for all real matrices A{t), necessarily 
rc < TT in (68), and more recently Moan and Niesen [175] have been able to 
prove that indeed Tc = tt if only real matrices are involved. 

In any case, it is important to remark that statement (D) is locally valid, but 
cannot be used to compute ft in the large. However, as we have seen, the other 
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statements need not depend on the validity of (D). In particular, if (B) and 
(C) are globally valid, one can still investigate many of the properties of fl 
even though one cannot compute it with the aid of (D). 



2.7.3 An improved radius of convergence 

The previous results on the convergence of the Magnus series have been estab- 
lished for n X n real matrices: if A{t) is a real nx n matrix, then (67) gives a 
condition for Y{t) to have a real logarithm. In fact, under the same condition, 
the Magnus scries (40) converges precisely to this logarithm, i.e., its sum Q(t) 
satisfies exp(0(t)) = Y{t) [175]. 

One should have in mind, however, that the original expansion was conceived 
by requiring only that A{t) be a linear operator depending on a real variable t 
in an associative ring (Theorem 1). The idea was to define, in terms of A, an 
operator Q{t) such that the solution of the initial value problem Y' = A(t)Y, 
Y{0) = I, for a second operator Y is given as F = expfi. The proposed 
expression for Q is an infinite series satisfying the condition that "its partial 
sums become Hermitian after multiplication by i if iA is a Hermitian operator" 
[158]. As this quotation illustrates, Magnus expansion was first derived in the 
context of quantum mechanics, and so one typically assumes that it is also 
valid when A{t) is a linear operator in a Hilbert space. Therefore, it might 
be desirable to have conditions for the convergence of the Magnus series in 
this more general setting. In [43] , by applying standard techniques of complex 
analysis and some elementary properties of the unit sphere, the bound Tc — tt 
has been shown to be also valid for any bounded normal operator A{t) in a 
Hilbert space of arbitrary dimension. Next we review the main issues involved 
and refer the reader to [43] for a more detailed treatment. 

Let us assume that A{t) is a bounded operator in a Hilbert space Ti, with 
2 < dim Ti < oo. Then we introduce a new parameter £ e C and denote by 
Y{t]e) the solution of the initial value problem 

^^^sAm yio) = i, (70) 

where now / denotes the identity operator in Tl. It is known that Y{t; e) is an 
analytic function of s for a fixed value of t. Let us introduce the set (Z C 
characterized by the real parameter 7, 



B^ = {eeC : \e\ [' \\A{s)\\ds < 'jj . 
Jo 



Here |{.|{ stands for the norm defined by the inner product on Ti, i.e., the 
2-norm introduced in subsection 1.2. 
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If t is fixed, the operator function ip{e) — logY{t; e) is well defined in when 
7 is small enough, say 7 < log 2, as an analytic function of e. As a matter of 
fact, this is a direct consequence of the results collected in section 2.7.2: if, 
in particular, /q ||A(s)||(is < log 2, the Magnus series corresponding to (70) 
converges and its sum Q{t;e) satisfies exp{Q{t; e)) — Y{t;e). In other words, 
the power series D,{t;e) coincides with (p(e) when \e\ Jq \\A{s)\\ds < log 2, and 
so the Magnus series is the power series expansion of ip{e) around £ = 0. 

Theorem 6 The function (p{e) = \ogY{t;s) is an analytic function of e in 
the set Bt^, with 

B^ = {eeC : |£| [' \\A(s)\\ds <7r}. 
Jo 

IfTi is infinite dimensional, the statement holds true ifV is a normal operator. 

In other words, 7 = vr. The proof of this theorem is based on some elementary 
properties of the unit sphere in a Hilbert space. Let us define the angle 
between any two vectors x ^ 0, y in H, Angja;, y} = a, < a < n, from 

Re{x,y) 
cos a — -p — — 1|- , 

where (•, •) is the inner product on H. This angle is a metric in S^, i.e., the 
triangle inequality holds in S^. 

The first property we need is given by the next lemma [174]. 

Lemma 7 For any x eH, x y^O, Ang{Y{t;e)x,x} < \e\ \\A{s)\\ds. 

Observe that if F is a normal operator in H, i.e., YY^ — Y^Y (in particular, if 
Y is unitary), then = for ell x E H and therefore Ang{Y^x, x} — 

Ang{Yx, x}. 

The second required property provides useful information on the location of 
the eigenvalues of a given bounded operator in H [171]. 

Lemma 8 Let T be a (bounded) operator on Ti. If AngjTx, x} < 7 and 
Ang{T'''a:, x} < 7 for any x 0, x E H, where denotes the adjoint operator 
ofT, then the spectrum ofT, o-{T), is contained in the set 

= 1^1 e^'^ e C : < 7}. 

Proof, (of Theorem 6). Let us introduce the operator T = Y{t; e), with e G 
B^, 7 < TT. Then by Lemma 7, Ang{Tx,x} < 7 for all x ^ 0, and thus, by 
Lemma 8, 

a(T) c A,. (71) 
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If dim Ti. — oo and we assume that Y{t; e) is a normal operator, then (71) also 
holds. 

Prom equation (70) in integral form, 

Y(t;e)^I + e C A{s)Yds, 
Jo 

one gets < 1 + |£| /q ||>l(s)|| II^IM*, and application of Gronwall's lemma 
[97] leads to 

\\Y{t;e)\\<exp(\e\ /j|A(.)|M.) . 
An analogous reasoning for the inverse operator also proves that 

\\Y-\t;e)\\<ew(^\e\J^\\A{s)\\dsy 

In consequence, 

||T||<e^ and < e^' . 

If A ^ G a{T), then |A| < ||T|| [108] and therefore |A| < e^. In addition, 
i e (7(T-^), so that |A| > e-T. Equivalently, 

a(T) C e C : e-T < \z\ < c^'} = G^. (72) 

Putting together (71) and (72), one has 

a(T) c n = A^. 

Now choose any value 7' such that 7 < 7' < tt (e.g., 7' = (7 + 7r)/2) and 
consider the closed curve P = dAy. Notice that the curve P encloses o"(T) in 
its interior, so that it is possible to define [74] the function (p{e) = logY{t;s) 
by the equation 

vie) = ^ / ^og z{zl - Y{t; e))-^ dz, (73) 

where the integration along P is performed in the counterclockwise direction. 
As is well known, (73) defines an analytic function of s in By [74] and the 
result of the theorem follows. ■ 

Theorem 9 Let us consider the differential equation Y' = A{t)Y defined in 
a Hilbert space 7i with Y{0) = I, and let A(t) be a bounded operator in Ti. 
Then, the Magnus series Q(t) = Y^'k'=i^kii) , with flk given by (51) converges 
in the interval i e [0, T) such that 

r\\A{s)\\ds<7: 
Jo 
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and the sum Q(t) satisfies expQ(t) = ^(0- statement also holds when 
Ti is infinite- dimensional if Y is a normal operator (in particular, if Y is 
unitary). 

Proof. Theorem 6 shows that \ogY{t; e) = ip{e) is a well defined and analytic 
function of s for 

|£| f \\Ais)\\ds < TT. 
^0 

It has also been shown that the Magnus series e) = J2kLi ^''^k{t): with 
given by (51), is absolutely converg cnt when \e\ £ \\A{s)\\ds < C = 1.0868... 
and its sum satisfies expQ{t;e) = Y{t;e). Hence, the Magnus series is the 
power series of the analytic function (p{e) in the disk \e\ < ^/ fQ\\A[s)\\ds. 
But ip{e) is analytic in D 5^ and the power series has to be unique. In 
consequence, the power series of ip{e) in has to be the same as the power 
series of ip{e) in S^, which is precisely the Magnus series. Finally, by taking 
£ = 1 we get the desired result. ■ 



2.7.4 Further discussion 



Theorem 9 provides sufficient conditions for the convergence of the Magnus 
series based on an estimate by the norm of the operator A. In particular, it 
guarantees that the operator Vt{t) in Y{t) — expQ(i) can safely be obtained 
with the convergent series ^k>i ^k{t) for < t < T when the terms ^^(t) are 
computed with (51). A natural question at this stage is what is the optimal 
convergence domain. In other words, is the bound estimate Tc = tt given by 
Theorem 9 sharp or is there still room for improvement? In order to clarify 
this issue, we next analyze two simple examples involving 2x2 matrices. 



Example 1. Moan and Niesen [175] consider the coefficient matrix 

A{t) 



(2 t\ 



(74) 



If we introduce, as before, the complex parameter e in the problem, the cor- 
responding exact solution Y{t;e) of (70) is given by 



Y{t;e)^ 



( „2£t _l_ -let 




-et 



1— £t 



(75) 



and therefore 



logF(t;£) 



2t g{t;e) 
-t 



, with g{t;e) 



J 



3(1 - e3^*) ' 
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The Magnus series can be obtained by computing the Taylor expansion of 
logY{t;e) around e = 0. Notice that the function g has a singularity when 
et = ^i, and thus, by taking s = 1, the Magnus series only converges up to t = 
gTT. On the other hand, condition /q ||A(s)||(is < TT leads to T ^ 1.43205 < fyr. 
In consequence, the actual convergence domain of the Magnus series is larger 
than the estimate provided by Theorem 9. □ 



Example 2. Let us introduce the matrices 



0-3, 




(76) 



and define 



A{t) 



(3X2 0<t<l 
aXi t>l 

with a, j3 complex constants. Then, the solution of equation (30) at t = 2 is 



Y{2) = e"-^i e^-^^ 



so that 



logy (2) = log(e"^i e^^2) = aX^ + 



2a(5 



^2, 



(77) 



1 - e-2" 

an analytic function if |a| < tt with first singularities at a = ±i7r. Therefore, 
the Magnus series cannot converge at i = 2 if |q;| > tt, independently of /? 7^ 0, 
even when it is possible in this case to get a closed-form expression for the 
general term. Specifically, a straightforward computation with the recurrence 
(49)- (50) shows that 



^ 1^,(2) = «Xi + /3X2 + ^(-l)^ 

n=l n=2 



-1 ^ -O5 



(n-1)! 



13X2. (78) 



If we take the spectral norm, then = 1 1-^^211 = 1 and 

'~^\\A{s)\\ds^\a\ + \(3\, 





so that the convergence domain provided by Theorem 9 is |q;| + < tt for 
this example. Notice that in the limit \j3\ — > this domain is optimal. □ 

From the analysis of Examples 1 and 2 we can conclude the following. First, 
the convergence domain of the Magnus series provided by Theorem 9 is the 
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best result one can get for a generic bounded operator A{t) in a Hilbert space, 
in the sense that one may consider specific matrices A{t), as in Example 2, 
where the series diverges for any time t such that /q ||yl(s)||(is > tt. Second, 
there are also situations (as in Example 1) where the bound estimate r — tt 
is still rather conservative: the Magnus series converges indeed for a larger 
time interval than that given by Theorem 9. This is particularly evident if one 
takes A{t) as a diagonal matrix: then, the exact solution Y{t\ e) of (70) is a 
diagonal matrix whose elements are non- vanishing entire functions of £, and 
obviously logF(t; e) is also an entire function of e. In such circumstances, the 
convergence domain \e\ fl \\A{s)\\ds < n for the Magnus series does not make 
much sense. Thus a natural question arises: is it possible to obtain a more 
precise criterion of convergence? In trying to answer this question, in [43] 
an alternative characterization of the convergence has been developed which 
is valid for n x n complex matrices. More precisely, a connection has been 
estabhshed between the convergence of the Magnus series and the existence 
of multiple eigenvalues of the fundamental matrix Y{t; e) for a fixed t, which 
we denote by Yt{e). By using the theory of analytic matrix functions, and in 
particular, of the logarithm of an analytic matrix function (such as is done e.g. 
in [242]), the following result has been proved in [43]: if the analytic matrix 
function Yt{e) has an eigenvalue po{eo) of multiplicity / > 1 for a certain Eq 
such that: (a) there is a curve in the e-plane joining e = with e = eo, and 
(b) the number of equal terms in \ogpi{eo), \ogp2{eo), ■ ■ ■ , logpi{eo) such that 
Pfc(^o) = po; ^ = 1, • • • ) ^ is less than the maximum dimension of the elementary 
Jordan block corresponding to po, then the radius of convergence of the series 
^ti^) = T,k>i^''^t,k verifying expQt(£) = Yt{e) is precisely r = |£o|- An 
analysis along the same line has been carried out in [230]. 

When this criterion is applied to Example 1, it gives as the radius of conver- 
gence of the Magnus series corresponding to equation (70) for a fixed t, 

oo 

nt{£) = T.^'^t,k, (79) 

k=l 

the value 

- I I - — 

To get the actual convergence domain of the usual Magnus expansion 
to take s — 1, and so, from (80), we get 27i/ (3t) = 1, or equivalently t 
i.e., the result achieved from the analysis of the exact solution. 

With respect to Example 2, one gets [43] 

vr 

" jaK^-l)' 

If we now fix £ = 1, the actual i-domain of convergence of the Magnus series 
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is 



Observe that, when t = 2, we get \a\ = vr and thus the previous result is 
recovered: the Magnus series converges only for \a\ < tt. 

It should also be mentioned that the case of a diagonal matrix A{t) is com- 
patible with this alternative characterization [43] . 



2.8 Magnus expansion and the BCH formula 

The Magnus expansion can also be used to get explicitly the terms of the 
series Z in 

Z = log(e^i e^'), 

Xi and X2 being two non commuting indeterminate variables. As it is well 
known [197], 

00 

Z = X,+X2 + J2 GniXi, X2), (81) 

n=2 

where Gn{Xi, X2) is a homogeneous Lie polynomial in Xi and X2 of grade 
n; in other words, Gn can be expressed in terms of Xi and X2 by addition, 
multiplication by rational numbers and nested commutators. This result is 
often known as the Bakcr-Campbell-Hausdorff (BCH) theorem and proves to 
be very useful in various fields of mathematics (theory of linear differential 
equations [158], Lie group theory [94], numerical analysis [99]) and theoretical 
physics (perturbation theory, transformation theory. Quantum Mechanics and 
Statistical Mechanics [142,238,239]). In particular, in the theory of Lie groups, 
with this theorem one can explicitly write the operation of multiplication in 
a Lie group in canonical coordinates in terms of the Lie bracket operation in 
its algebra and also prove the existence of a local Lie group with a given Lie 
algebra [94]. 



If Xi and X2 are matrices and one considers the piecewise constant matrix- 
valued function 



X2 



0<t<l 
l<t<2 



(82) 



then the exact solution of (30) at t = 2 is Y(2) = e ^ e 2. By computing 
y(2) = e^(^) with recursion (51) one gets for the first terms 
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^^i(2) 
1^2(2) 



^[X,,[X,,X,]]- ^[X,,[X,,X,]] 
i-[Xi,[X2,[X2,Xi]]] 



(83) 



^3(2) 



1^4(2) 



In general, each G„(Xi, X2) is a linear combination of the commutators of the 



form [Vi, [V2, [Vn-i, Vn] . . .]] with Vi E {Xi, X2} for 1 < i < n, the coeffi- 



cients being universal rational constants. This is perhaps one of the reasons 
why the Magnus expansion is often referred to in the literature as the contin- 
uous analogue of the BCH formula. As a matter of fact, Magnus proposed a 
different method for obtaining the first terms in the series (40) based on (81) 



Now we can apply Theorem 9 and obtain the following sharp bound. 

Theorem 10 The Baker-Camphell-Hausdorff series in the form (81) con- 
verges absolutely when \\Xi\\ + ||X2|| < tt. 

This result can be generalized, of course, to any number of non commuting 
operators Xi, X2, . . . , X^. Specifically, the series 



converges absolutely if -h IIX2II -!-■■■ + < tt. 
2.9 Preliminary linear transformations 

To improve the accuracy and the bounds on the convergence domain of the 
Magnus series for a given problem, it is quite common to consider first a linear 
transformation on the system in such a way that the resulting differential 
equation can be more easily handled in a certain sense to be specified for each 
problem. To illustrate the procedure, let us consider a simple example. 

Example. Suppose we have the 2x2 matrix 



where Xi and X2 are given by (76) and a and (5 are complex functions of 
time, a, /3 : R — > C Then the exact solution of Y' = A{t)Y , y (0) = 7 is 



[158]. 



Z = log(e^i e^2 ... gX,^^ 



A{t) = a{t)Xi + (3{t)X, 



2, 



(84) 




(85) 
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Let us factorize the solution as Y(t) — Yo{t)Yi{t), with Yo{t) the solution of 
the initial value problem defined by 

(a(t) \ 

yj = Ao{t) Yo Ao{t) = a{t)Xi = (86) 

and Yq{0) = I. Then, the equation satisfied by Yi is 

^ e^./'o"^^)'^^^ 



y/ = Ai {t)Yi , with Ai = Y^-^A Yo - Aq ^ 



V / 



(87) 

so that the first term of the Magnus expansion applied to (87) already provides 
the exact solution (85). □ 

This, of course, is not the typical behaviour, but in any case if the transfor- 
mation Yq in the factorization Y{t) = Yo{t)Yi{t) is chosen appropriately, the 
first few terms in the Magnus series applied to the equation satisfied by of Yi 
give usually very accurate approximations. 

Since this kind of preliminary transformation is frequently used in Quantum 
Mechanics, we speciahze the treatment to this particular setting here and 
consider equation (4) instead. In other words, we write (30) in the more con- 
ventional form of the time dependent Schodinger equation 

^^mUit), (88) 

where H = H/{ih), h is the reduced Planck constant, H is the Hamiltonian 
and U corresponds to the evolution operator. 

As in the example, suppose that H can be split into two pieces, H — Hq+eHi, 
with Hq a solvable Hamiltonian and e ^ 1 a small perturbation parameter. In 
such a situation one tries to integrate out the Hq piece so as to circumscribe the 
approximation to the Hi piece. In the case of equation (88) this is carried out 
by means of a linear time-dependent transformation. In Quantum Mechanics 
this preliminary linear transformation corresponds to a new evolution picture, 
such as the interaction or the adiabatic picture. 

Among other possibilities we may factorize the time-evolution operator as 

U{t)^Git)UG{t)G^{0), (89) 

where G{t) is a linear transformation whose purpose is to be defined yet. In 
the new G-Picture, the corresponding time-evolution operator Uq obeys the 



38 



equation 

U'ait) = HcmGit), Hcit) = G\t)H{t)G{t) - G\t)G\t). (90) 

The choice of G depends on the nature of the problem at hand. There is 
no generic formal recipe to find out the most appropriate G. In the spirit 
of canonical transformations of Classical Mechanics, one should built up the 
very Uq perturbatively. However, the aim here is different because G is defined 
from the beginning. Two rather common choices are: 

• Interaction Picture. It is well suited when i?o(^) is diagonal in some basis, 
or else, it is constant. In that case 



G{t) = exp ^o(r)dr) (91) 



so that 



Hcit) = e exp (- Ho{,r)dT^ Hi{t) exp Ho{r)dT^ . (92) 

• Adiabatic Picture. A time scale of the system much smaller than that of 
the interaction defines an adiabatic regime. For instance, suppose that the 
Hamiltonian operator H(t) depends smoothly on time through the variable 
r = t/T, where T determines the time scale and T ^ oo. Then the quantum 
mechanical evolution of the system is described by dU/dt = H{st)U, with 
e = 1/T <^ 1, or equivalently 

with T = et. In this case the appropriate transformation is a G that renders 
H{t) instantaneously diagonal, i.e., 

G\t)H{t)G{t)^E{t)^dieig[E,{t),E2{t),...]. (94) 

The term G^G" of the new Hamiltonian in (90) is, under adiabatic condi- 
tions, very small. Its main diagonal generates the so-called Berry, or geo- 
metric, phase [14]. 

Both types of G do not exclude mutually, but they may be used in succes- 
sion. As a matter of fact, corrections to the adiabatic approximation must 
be followed by the former one. In turn, an adiabatic transformation may be 
iterated, as proposed by Garrido [92] and Berry [15] . 

In section 4 we shall use extensively these preliminary linear transformations 
on several standard problems of Quantum Mechanics to illustrate the practical 
features of the Magnus expansion. 



39 



2.10 Exponential product representations 



In contrast to Magnus expansion, much less attention has been paid to so- 
lutions of (30) in the form of a product of exponential operators. Both ap- 
proaches are by no means equivalent, since in general the operators Q„ do not 
commute with each other. For instance, for a quantum system as in equation 
(88), the ansatz U = nexp($„) (where are skew-Hermitian operators to 
be determined) is an alternative to the Magnus expansion, also preserving the 
unitarity of the time-evolution operator. One such procedure was devised by 
Fer in 1958 in a paper devoted to the study of systems of linear differential 
equations [85] . Although the original result obtained by Fer was cited and ex- 
plicitly stated by Bellman [13, p. 204], sometimes it has been misquoted as a 
reference for the Magnus expansion [11]. On the other hand, Wilcox associated 
Fer's name with an interesting alternative infinite product expansion which 
is indeed a continuous analogue of the Zassenhaus formula [239] (something 
also attributed to the Fer factorization [159, p. 372]). This however also led 
to some confusion since his approach is in the spirit of perturbation theory, 
whereas Fer's original one was essentially nonperturbative. The situation was 
clarified in [136], where also some applications to Quantum Mechanics were 
carried out for the first time. 

In this section wc discuss briefly the main features of the Fer and Wilcox 
expansions, and how the latter can be derived from the successive terms of 
the Magnus series. This will make clear the different character of the two ex- 
pansions. We also include some details on the factorization of the solution 
proposed by Wei and Norman [236,237]. Finally we provide another inter- 
pretation of the Magnus expansion as the continuous analogue of the BCH 
formula in linear control theory. 



2.10.1 Fer method 

An intuitive way to introduce Fer formalism is the following [119]. Given the 
matrix linear system Y' = A{t)Y, Y{0) — I, we know that 

Y{t)=exp{F,{t)) (95) 

is the exact solution if A commutes with its time integral Fi{t) — jQA{s)ds, 
and Y(t) evolves in the Lie group QUA lies in its corresponding Lie algebra 
Q. If the goal is to respect the Lie-group structure in the general case, we need 
to 'correct' (95) without loosing this important feature. 

Two possible remedies arise in a quite natural way. The first is just to seek a 
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correction A{t) evolving in the Lie algebra g so that 

=exp (Fi(t) + A(t)). 

This is nothing but the Magnus expansion. Alternatively, one may correct 
with Yi{t) in the Lie group Q, 

Y{t)=exp{F,{t))Y,{t). (96) 

This is precisely the approach pursued by Fer, i.e. representing the solution 
of (30) in the factorized form (96), where (hopefully) Yi will be closer to the 
identity matrix than Y at least for small t. 

The question now is to find the differential equation satisfied by Yi. Substi- 
tuting (96) into equation (30) we have 

so that, taking into account the expression for the derivative of the exponential 
map (Lemma 2), one arrives easily at 

y; = A,{t)Y, Y,(0) = I, (98) 

where 

(t) = e-^i A e^i - dx e'^-^i A e^^^ . (99) 
Jo 

The above procedure can be repeated to yield a sequence of iterated matrices 
Ak- After n steps we have the following recursive scheme, known as the Fer 
expansion: 

Y = e^' c^' ■ ■ ■ c^" Yr, (100) 
y^ = ^„(t)F„ Y^{0) = I, J = 1,2,... 

with Fn{t) and A„(f) given by 



Fn+i{t) = /* Ar,{s)ds Ao{t) = A{t), n = 0, 1, 2... 
Jo 

= e-^^+' An e^"+i - f dx e'^-^^+i ^„ e^-^"+i 
Jo 

= [ dx r du e-(^-")^"+i e(^-")^"+i (101) 

JO Jo 



Ejj^_^dk,Mn), 0,1,2.. 
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When after n steps we impose Yn — I we are left with an approximation to 
the exact solution Y{t). 

Inspection of the expression of An+i in (101) reveals an interesting feature of 

the Fcr expansion. If wc assume that a perturbation parameter e is introduced 
in A, i.e. if wc substitute A by sA in the formalism, since is of the same 
order in e as An, then an elementary recursion shows that the matrix An starts 
with a term of order e^" (correspondingly the operator Fn contains terms of 
order and higher). This should greatly enhance the rate of convergence 
of the product in equation (100) to the exact solution. 

It is possible to derive a bound on the convergence domain in time of the 
expansion [21]. The idea is just to look for conditions on A{t) which insure 
F„ — > as n — > oo. As in the case of the Magnus expansion, we take A{t) to be 
a bounded matrix with ||A(i)|| < k{t) = ko{t). Fer's algorithm, equations (100) 
and (101), provides then a recursive relation among corresponding bounds 
kn{t) for ||74„(t)||. If we denote Kn{t) = /q kn{s)ds, we can write this relation 
in the generic form kn+i = f{kn, Kn), which after integration gives 

K„+,^M{K„). (102) 



The question now is: When is Kn — as n — oo? This is certainly so if is 
a stable fixed point for the iteration of the mapping M and Kq is within its 
basin of attraction. To see when this is the case we have to solve the equation 
^ = M{^) to find where the next fixed point hes. Let us do it explicitly. By 
taking norms in the recursive scheme (101) we have 



ll^n+ill < f'dx r du e2(i-»)^"||K,F„ 
Jo Jo 



+1JII > 

which can be written as ||A„+i|| < A;„+i,with 

_ l-e=^^"(l-2i^^) dK^ 
~ 2Kn dt 

and consequently Kn+i is given by eq. (102) with 

M(A„) = f dx. (103) 

That is the mapping we have to iterate. It is clear that ^ = is a stable 
fixed point of M. The next, unstable, fixed point is ,^ = 0.8604065. So we can 
conclude that we have a convergent Per expansion at least for values of time 
t such that ^ 

/ \\A{s)\\ds < Ko{t) < 0M0m5. (104) 
Jo 

Notice that the bound for the convergence domain provided by this result is 
smaller than the corresponding to the Magnus expansion (Theorem 9). 
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2.10.2 Wilcox method 



A more tractable form of infinite product expansion lias been devised by 
Wilcox [239] in analogy with the Magnus approach. The idea, as usual, is 
to treat e in 

y = eA{t)Y, Y{0) = I (105) 

as an expansion parameter and to determine the successive factors in the 
product 

Y{t) = e^i e^2 e^3 . . . (io6) 
by assuming that Wn is exactly of order s". Hence, it is clear from the very 
beginning that the methods of Fer and Wilcox give rise indeed to completely 
different infinite product representations of the solution Y{t). 

The explicit expressions of Wi, W2 and are given in [239]. It is noteworthy 
that the operators Wn can be expressed in terms of Magnus operators ftk, for 
which compact formulae and recursive procedures are available. To this end we 
simply use the Baker-Campbell-Hausdorff formula to extract formally from 
the identity 

e^^i . . . = ^n,+n,+ns+- ^ (io7) 

terms of the same order in e. After a straightforward calculation one finds for 
the first terms 



w^i = fii, W2 = n2, Ws = ns-^[ni,n2i (los) 

W^^n^- li^u^s] + li^u [^1:^2]]: etc. (109) 



The main interest of the Wilcox formalism stems from the fact that it provides 
explicit expressions for the successive approximations to a solution represented 
as an infinite product of exponential operators. This offers a useful alternative 
to the Fer expansion whenever the computation of Fn from equation (101) is 
too cumbersome. We note in passing that to first order the three expansions 
yield the same result: Fi = l^i = Qi. 



2.10.3 Wei-Norman factorization 

Suppose now that A and Y in equation (30) are linear operators and that A{t) 
can be expressed in the form 

m 

A{t) = Y,Ui{t)Xi, m finite, (110) 

i=l 

where the Ui{t) are scalar functions of time, and Xi,X2, are time- 

independent operators. Furthermore, suppose that the Lie algebra Q generated 
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by the Xi is of finite dimension / (this is obviously true if A and Y are finite 
dimensional matrix operators). Under these conditions, if Xi,X2, . . . ,Xi is a 
basis for q, the Magnus expansion allows to express the solution locally in the 
form Y{t) = exp(X]'=i /i(t)^i)- Wei and Norman, on the other hand, show- 
that there exists a neighborhood oit — in which the solution can be written 
as a product [236,237] 

Y{t) = eMgiiW eM92{t)X2) ■ ■ -eMMXi), (111) 

where the gi{t) are scalar functions of time. Moreover, the gi{t) satisfy a set of 
nonlinear differential equations which depend only on the Lie algebra g and 
the Ui{tys. These authors also study the conditions under which the solution 
converges globally, that is, for all t. In particular, this happens for all solvable 
Lie algebras in a suitable basis and for any real 2x2 system of equations [237]. 

In the terminology of Lie algebras and Lie groups, the representation Y{t) — 
^^p{J2i=i fi{t)Xi) corresponds to the canonical coordinates of the first kind, 
whereas equation (111) defines a system of canonical coordinates of the second 
kind [12,94,197]. 

This class of factorization has been used in combination with the Fer expan- 
sion to obtain closed-form solutions of the Cauchy problem defined by certain 
classes of parabolic linear partial differential equations [41]. When the algo- 
rithm is applied, the solution is written as a finite product of exponentials 
depending on certain ordering functions for which convergent approximations 
are constructed in explicit form. 

Notice that the representation (111) is clearly useful when the spectral prop- 
erties of the individual operators Xi are readily available. Since the Xi are 
constant and often have simple physical interpretation, the evaluation of the 
eigenvalues and eigenvectors can be done once for all times, and this may facil- 
itate the computation of the exponentials. This situation arises, in particular, 
in control theory [12]. The functions Ui{t) are known as the controls, and the 
operator Y(t) acts on the states of the system, describing how the states are 
transformed along time. 



2.11 The continuous BCH formula 

When applied to the equation Y' = A{t)Y with the matrix A(t) given by 
(110), the Magnus expansion adopts a particularly simple form. Furthermore, 
by making use of the structure constants of the Lie algebra, it is relatively 
easy to get explicit expressions for the canonical coordinates of the first kind 
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fi{t). Let us illustrate the procedure by considering the particular case 

A{t) = Ui{t)Xi + U2{t)X2. 

Denoting by ai{t) — jQUi{s)ds, and, for a given function //, 

a straightforward calculation shows that the first terms of Q in the Magnus 
expansion can be written as 

Q{t)=f3,{t)X, + l32{t)X2 + A2(t)[Xi,X2] + Ai2(t)[Xi, [Xi,X2]] 

+/9212(t)[X2,[Xi,X2]] + --- (112) 

where 



/3i2 = ^(/«2-^«i), (113) 

^T2^l/^-L "^"^^ - i^/i i - /i / 

/3212 = ^(^ aia2 - / "2) + ^(/J «2 - ai). 

Taking into account the structure constants of the particular finite dimensional 
Lie algebra under consideration, from (112) one easily gets the functions fi{t). 
In the general case, (112) allows us to express f2 as a linear combination of 
elements of a basis of the free Lie algebra generated by Xi and X2. In this 
case, the recurrence (49)- (50) defining the Magnus expansion can be carried 
out only with the nested integrals 

(i) = (^^ • • • ^ 1) {t) ^ j'^ j'j ■■■ J^' J^' {t,) ■■■Ui, {h)dh ■■■dts 

(114) 

involving the functions Ui{t) and U2{t). Thus, for instance, the coefficients in 
(113) can be written (after successive integration by parts) as 



/3j Q; j , i 1)2, 

/3i2 = ^(^q;2 -^Qii) = ^(q;2i - q;i2), 



/3ii2 = J(/ 1 0L1+ I I a2-2 I f ai) = ^(aii2 + "211 - 2ai2i), 



6 J2 Jl Jl Jl Jl J2 6 

i,^ r r r r r r . i 



^212 = -{2 f fa2-f / tti - / f 02) = -(2q;212 - ai22 - Q;22i)- 
J2J1 J2J2 J1J2 
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The series (112) expressed in terms of the integrals (114) is usually referred to 
as the continuous Bakcr-Campbell-HausdorfT formula [132,186] for the linear 
case. We will generalize this formalism to the nonlinear case in the next section. 



3 Generalizations of the Magnus expansion 

In view of the attractive properties of the Magnus expansion as a tool to 
construct approximate solutions of non-autonomous systems of linear ordinary 
differential equations, it is hardly surprising that several attempts have been 
made along the years either to extend the procedure to a more general setting 
or to manipulate the series to achieve further improvements. In this section we 
review some of these generalizations, with special emphasis on the treatment 
of nonlinear differential equations. 

First we reconsider an iterative method originally devised by Voslamber [232] 
for computing approximations Q^'^\t) in Y{t) — exp{Q{t)) for the linear equa- 
tion y = A{t)Y. The resulting approximation may be interpreted 
summation of terms in the Magnus series and possesses interesting features 
not shared by the corresponding truncation of the conventional Magnus ex- 
pansion. Then we adapt the Magnus expansion to the physically relevant case 
of a periodic matrix A{t) with period T which incorporates in a natural way 
the structure of the solution ensured by the Floquet theorem. Next we go 
one step further and generalize the Magnus expansion to the so-called non- 
linear Lie equation Y' = A{t, Y)Y. Finally, we show how the procedure can 
be applied to any nonlinear explicitly time-dependent differential equation. 
Although the treatment is largely formal, in section 5 we will see that it is 
of paramount importance for designing new and highly efficient numerical in- 
tegration schemes for this class of differential equations. We particularize the 
treatment to the important case of Hamiltonian systems and also establish 
an interesting connection with the Chen-Fliess series for nonlinear differential 
equations. 



3.1 Voslamber iterative method 

Let us consider equation (30) when there is a perturbation parameter e in the 
(in general, complex) matrix A, i.e., equation (105). Theorem 9 guarantees 
that, for sufficiently small values of t, Y{t;e) — expfl{t;e), where 

oo 

n{t;e) = Y. sennit). (115) 

n=l 
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The advantages of this representation and the approximations obtained when 
the series is truncated have been sufficiently recognized in the treatment done 
in previous sections. There is, however, a property of the exact solution not 
shared by any truncation of the series (115) which could be relevant in certain 
physical applications: {1 / e) Y^^=i £"^n{t) with m > 1 is unbounded for e 
oo even when n(t,£)/e is bounded uniformly with respect to e under rather 
general assumptions on the matrix A{t) [232]. Notice that this is the case, in 
particular, for the adiabatic problem (93). 

When Schur's unitary triangularization theorem [107] is applied to the exact 
solution Y{t;e) one has 

T, = U^Y U, (116) 

where is an upper triangular matrix and U is unitary. In other words, Y 
is unitarily equivalent to an upper triangular matrix T^. Differentiating (116) 
and using (105) one arrives at 

= sU^'AUT, + [t„ U^U' 

Since the second term on the right hand side is not upper triangular, it follows 
at once that 



where the subscript A denotes the upper triangular part (including terms on 
the main diagonal) of the corresponding matrix. Taking into account (116) 
one gets 

n{t, e) = eU (^J^ {U^AU)^ds^ (117) 

Considering now the Frobenius norm (which is unitarily invariant, section 1.2) 
of both sides of this equation, one has 



|J1||f= kl 



f\u^AU)^ds <\e\ f' \\{U''AU)4Fds 

Jo F Jo 



< \e\ 



WU^AUyds 



\\A\\Fds, 



(118) 



If the the spectral norm is considered instead, from inequalities (28), (29) and 
(118), one concludes that 



\n\\2 < ^rank(A) |£| \\A\\2 ds. 



In any case, what is important to stress here is that for the exact solution 
fl{t;e)/6 is bounded uniformly with respect to the e parameter. Voslambcr 
proceeds by deriving an algorithm for generating successive approximations 
of Y{t; e) = exp{n{t; e)) which, contrarily to the direct series expansion (115), 
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preserve this property. His point of departure is to get a series expansion for 
the so-called dressed derivative of VL [193] 

r = e"/2 e-"/^ (119) 

This is accomplished by inserting (39) in (119). Specifically, one has 



sinhO/2^^'^^-„eo n\ ^^"^^"^^ 



and finally [232,193] 



oo Ol— n 1 

r=E ^B^ad^sA), (120) 

where, as usual, denote Bernoulli numbers. In order to express F as a power 
series of e one has to insert the Magnus series (115) into eq. (120). Then we 
get 

oo 

r{t;e) = Y.^''^n{t), (121) 

n=l 

where the terms can be expressed as a function of Jl^ with k < n — 2 
through the recursive procedure [193] 

Ti = A, F2 = 0, (122) 

n— 1 

^n = Yl 9 E adn,^ adf^,^ • • • adn, . A, n>3. 

j=2 ki + --- + kj=n-l 
ki>l,...,kj>l 

Here 

_ 2^-^ - 1 
Cj = - Bj, 

with C2j+i — 0, C2 — —1/24, C4 — 7/5760, etc. In particular. 



r3 = -^[fii,[fii,A]] 

r4 = -^([Oi,[Q2,A]] + [Q2,[^^l,A]]). 

Now, from the definition of F, eq. (119), we write 

n' = e-"/2 F e"/^ 
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which, after integration over t, can be used for constructing successive approx- 
imations to Q once the terms r„ are known in terms of 0^, k < n — 2. Thus, 
the nth approximant ^2*^"^ is defined by 

QH(i) = /*e-|f^^"-^n^)rW(s) ei"^""'(^)ds, n= 1,2,... (123) 
Jo 

where the e dependence has been omitted by simphcity and F^") = Sfe=i ^'^^k, 
— O. The first two approximants read exphcitly 

n(^\t,6)=6Qi{t) =e t A{s)ds 

Jo 

/*e-5^'''(^'-)>l(s) e5^''H^'-)ds. (124) 
Jo 

In this approach the solution is approximated by Y{t) ~ exp(r2*^")). Observe 
that contains contributions from an infinity of orders in s, whereas the 
nth term in the Magnus series (115) is proportional to Furthermore, fi^") 
contains Y.k=i^''^k c-^d also higher powers e'^ (m > n). In particular, one 
gets easily 

From the structure of the expression (123) it is also possible to find the asymp- 
totic behaviour of ^2*^"^^ ^ 3) for e ^ oo and prove that it remains bounded 
[232], just as the exact solution does. This property of the Voslamber iterative 
algorithm may lead to better approximations of Y{t) when the parameter e is 
not very small, since in that case Q^'^^s is expected to remain close to Jl/s, 
as shown in [193]. 

3.2 Floquet-Magnus expansion 

We now turn our attention to a specific case of equation (30) with important 
physical and mathematical applications, namely when the (complex) n x n 
matrix-valued function A{t) is periodic with period T. Then further infor- 
mation is available on the structure of the exact solution as is given by the 
celebrated Floquet theorem, which ensures the factorization of the solution in 
a periodic part and a purely exponential factor. More specifically, 

Y(t) ^ P(t) ex.p{tF), (125) 

where F and P are n x n matrices, P{t) = P{t + T) for all t and F is con- 
stant. Thus, albeit a solution of (30) is not, in general, periodic the departure 
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from periodicity is determined by (125). This result, when apphed in quan- 
tum mechanics, is referred to as Bloch wave theory [29,89]. It is widely used 
in problems of solid state physics where space-periodic potentials are quite 
common. In Nuclear Magnetic Resonance this structure is exploited as far as 
either time-dependent periodic magnetic fields or sample spinning are involved 
[82]. Asymptotic stability of the solution Y(t) is dictated by the nature of the 
eigenvalues of F, the so-called characteristic exponents of the original periodic 
system [101]. 

An alternative manner of interpreting equation (125) is to consider the piece 
P{t), provided it is invertible, to perform a transformation of the solution in 
such a way that the coefficient matrix corresponding to the new representation 
has all its matrix entries given by constants. Thus the piece exp(tF) in (125) 
may be considered as an exact solution of the system (30) previously moved to 
a representation where the coefficient matrix is the constant matrix F [242]. 
The t-dependent change of representation is carried out by P{t). Connecting 
with section 2.9, P{t) is the appropriate preliminary linear transformation 
for periodic systems. Of course, Floquet theorem by itself gives no practical 
information about this procedure. It just states that such a representation 
does exist. In fact, a serious difficulty in the study of differential equations 
with periodic coefficients is that no general method to compute either the 
matrix P[t) or the eigenvalues of F is known. 

Mainly, two ways of exploiting the above structure of Y{t) are found in the 
literature [147]. The first one consists in performing a Fourier expansion of 
the formal solution leading to an infinite system of linear differential equations 
with constant coefficients. Thus, the t-dependent finite system is replaced with 
a constant one at the price of handling infinite dimension. Resolution of the 
truncated system furnishes an approximate solution. The second approach is 
of perturbative nature. It deals directly with the form (125) by expanding 

oo oo 

P{t)^Y.Pnit), F^Y.Fn. (126) 

n=l n=l 

Every term F„ in (126) is fixed so as to ensure Pn{t) = Pn(t + T), which in 
turn guarantees the Floquet structure (125) at any order of approximation. 

Although the Magnus expansion, such as it has been formulated in this work, 
does not provide explicitly the structure of the solution ensured by Floquet 

theorem, it can be adapted without special difficulty to cope also with this 
situation. The starting point is to introduce the Floquet form (125) into the 
differential equation Y' = A{t)Y. In that way the evolution equation for P is 
obtained: 

P'{t) = A{t)P{t) - P{t)F, P(0) = I. (127) 
The constant matrix F is also unknown and we will determine it so as to 
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ensure P(t + T) = P{t). Now we replace the usual perturbative scheme in 
equation (126) with the exponential ansatz 

P{t) = exp(A(t)), A(0) = O. (128) 

Obviously, A{t + T) — A{t) so as to preserve periodicity. Now equation (127) 
conveys 

^ cxp(A) = A exp(A) - exp(A)F, (129) 
from which, as with the conventional Magnus expansion, it follows readily that 

oo r) 

^' = Y.^^dHA + {-lf+'F). (130) 

fe=0 

This equation is now, in the Floquet context, the analogue of Magnus equation 
(39). Notice that if we put F = O then (39) is recovered. The next move is to 
consider the series expansions for A and F 



A{t) = J2Ak{t), F=Y.Fk, (131) 

k=l k=l 

with Ak{0) = O, for all k. Equating terms of the same order in (130) one gets 
the successive contributions to the series (131). Therefore, the explicit ansatz 
we are propounding reads 

Y{t) = exp l^f; Afc(t) j exp (^t f] F,. j . (132) 

This can be properly referred as the Floquet-Magnus expansion. 

Substituting the expansions of equation (131) into (130) and equating terms 
of the same order one can write 

n— 1 r) 

a; = E -f {^n\t) + (-l)^+^r(^)(^)) in > 1). (133) 

j=0 J- 

The terms Wj^^\t) may be obtained by a similar recurrence to that given in 
equation (49) 

Wi^^ = [Am, Wtr^] {l<J<n- 1), 

m=l 

(134) 

wi^'^^A, lyW^o (n>l). 
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whereas the terms T^^^ (t) obey to the recurrence relation 



n-j 



Tjf^ = E A.., T}f_-^^\ (1 < J < n - 1), 

m=l 

(135) 

(n > 0). 

Every F„ is fixed by the condition A„(t + T) = A„(t). An outstanding feature 
is that F„ can be determined independently of A„(i) as the solution Y{t) — 
P(i)exp(iF) shrinks to Y{T) — exp(T'F). Consequently, the conventional 
Magnus expansion Y(t) — exp(Q(t)) computed at t — T must furnish 

F„ = for all n. (136) 

The first contributions to the Floquet-Magnus expansion read exphcitly 



Ai(t)= / A(x) dx-tFi, 
Jo 

Fi = ^ [ A{x) dx, 

K2{t)^\ [\a{x) + F,,A^{x)] dx-tF2, (137) 
2 Jo 

F2 = ^ £[Aix) + F,,A^ix)] dx. 



Moreover, from the recurrence relations (134) and (135) it is possible to obtain 
a sufficient condition such that convergence of the series J2 is guaranteed 
in the whole interval t E [0, T] [45]. In fact, one can show that absolute con- 
vergence of the Floquet-Magnus series is ensured at least if 



/ |U(i) II dt < = 0.20925. (138) 
Jo 

Notice that convergence of the series J^Fn is already guaranteed by (136) and 
the discussion concerning the conventional Magnus expansion in subsections 
2.7.2 and 2.7.3. The bound in the periodic Floquet case turns out to be 
smaller than the corresponding bound Tc = tt in the conventional Magnus 
expansion. At first sight this could be understood as an impoverishment of the 
result. However it has to be recalled that, due precisely to Floquet theorem, 
once the condition is fulfilled in one period convergence is assured for any 
value of time. On the contrary, in the general Magnus case the bound gives 
always a running condition. 
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3.3 Magnus expansions for nonlinear matrix equations 

It is possible to extend the procedure leading to the Magnus expansion for 
the linear equation (30) and obtain approximate solutions for the nonlinear 
matrix equation 

Y' = A{t, Y)Y, F(0) = Yoeg, (139) 

where ^ is a matrix Lie group, A : x Q — > g and q denotes the correspond- 
ing Lie algebra. Equation (139) appears in relevant physical fields such as rigid 
body mechanics, in the calculation of Lyapunov exponents {Q = SO(n)) and 
other problems arising in Hamiltonian dynamics {Q = Sp(n)). In fact, it can 
be shown that every differential equation evolving on a matrix Lie group Q can 
be written in the form (139) [119]. Moreover, the analysis of generic differen- 
tial equations defined in homogeneous spaces can be reduced to the Lie-group 
equation (139) [184]. 

In [44] a general procedure for devising Magnus expansions for the nonlinear 
equation (139) is introduced. It is based on applying Picard's iteration on the 
associated differential equation in the Lie algebra and retaining in each itera- 
tion the terms necessary to increase the order while maintaining the explicit 
character of the expansion. The resulting methods are thus explicit by design 
and are expressed in terms of integrals. 

As usual, the starting point in the formalism is to represent the solution of 
(139) in the form 

Y(t)^eM^(t,Yo))Yo. (140) 
Then one obtains the differential equation satisfied by O,: 

n' = d expf^^ {A{t, e^ Fo)) , 1^(0) = O, (141) 

where dexp^^ is given by (38). Now, as in the linear case, one can apply 
Picard's iteration to equation (141), giving instead 

■'O fc=0 

The next step in getting explicit approximations is to truncate appropriately 

the (iexp~^ operator in the above expansion. Roughly speaking, when the 
whole series for dexp^^ is considered, the power series expansion of the iterate 
function VL^^^t)., k > 1, only reproduces the expansion of the solution il[t) up 
to certain order, say m. In consequence, the (infinite) power series of fl^''\t) 
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and differ in terms The idea is then to discard in all 

terms of order greater than m. This of course requires careful analysis of each 
term in the expansion. For instance, fi™ = O implies that (f^I^l)' = ^(t, Fq) 
and therefore 

= A(s,Yo)ds ^ n(t,Yo) + 0(t^). 
Jo 

Since 

A(s,e^'''(^)Fo) = A(0,lo) + 0(5) 

it follows at once that 

-I fp%),A{s,e'''''(^)Yo)]ds^O{t'). 
2 Jo 

When this second term in ^^[^'(t) is included and i^'^l is computed, it turns out 
that fi'^l reproduces correctly the expression of f)!^! up to 0{t'^). Therefore we 
truncate rfexp^^ at the k = term and take 

Qpi(i)= rA(s,e"'''(^)yo)rfs. 

Jo 

With greater generality, we let 

QW(t)= [' A(s,Yo)ds (142) 
Jo 

fc=o -'^ 

and take the approximation Vt{t) ~ ^['"'(t). This results in an explicit ap- 
proximate solution that involves a linear combination of multiple integrals of 
nested commutators, so that ^["'^(t) G for all m > 1. It can also be proved 
that once inserted in (140), provides an explicit approximation yt™'l(t) 

for the solution of (139) that is correct up to terms 0{t'^^^) [44]. In addition, 
qH(^) reproduces exactly the sum of the first m terms in the Vt series of the 
usual Magnus expansion for the linear equation Y' = A(t)Y. It makes sense, 
then, to regard the scheme (142) as an explicit Magnus expansion for the 
nonlinear equation (139). 

This procedure can be easily adapted to construct an exponential representa- 
tion of the solution for the differential system 

Y' = [Ait, Y),Y], y(0) = Fo e Sym(n). (143) 

Here Sym(n) stands for the set of n x n symmetric real matrices and the 
(sufficiently smooth) function A maps R+ x Sym(n) into so(n), the Lie alge- 
bra oi n X n real skew-symmetric matrices. It is well known that the solution 
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itself remains in Sym(n) for all t > 0. Furthermore, the eigenvalues of Y{t) 
are independent of time, i.e., Y{t) has the same eigenvalues as Yq. This re- 
markable qualitative feature of the system (143) is the reason why it is called 
an isospectral flow. Such flows have several interesting applications in physics 
and applied mathematics, from molecular dynamics to micromagnetics to lin- 
ear algebra [39]. 

Since Y{t) and Y{0) share the same spectrum, there exists a matrix function 
Q{t) e SO(n) such that Y{t)Q{t) = Q{t)Y{0) or, equivalently, 

Y{t)^Qit)YoQ'^{t). (144) 

Then, by inserting (144) into (143), it is clear that the time evolution of Q{t) 
is described by 

Q' = A{t, QYoQ^) Q, g(0) = /, (145) 

i.e., an equation of type (139). Yet there is another possibility: if we seek 
the orthogonal matrix solution of (145) as Q{t) — exp{il{t)) with fl skew- 
symmetric, 

Y{t) = e"W Yo e-"W, t > 0, fl{t) e 5o{n), (146) 
then the corresponding equation for ft reads 

n' = dexp^^ [A{t,e^Yoe-^)'), Q(0) = O. (147) 

In a similar way as for equation (141), we apply Picard's iteration to (147) and 
truncate the ciexp"^ series at /c = m — 2. Now we can also truncate consistently 
the operator 

Adnyo = e"yoe-" = e^"yo 

and the outcome still lies in so(n). By doing so, we replace the computation 
of one matrix exponential by several commutators. 

In the end, the scheme reads 

QW(t)= f' A(s,Yo)ds 
Jo 

m— 1 2 

Qm-iit) = E 7[a4[^-i](,)yo (148) 

m— 2 75 .1 

qH(^) = E it / ad^[_,](,3A(s, em-iis))ds, m > 2 
k=o •'^ 

and, as before, one has VL{t) = (]H(t) + ^(t^+i). Thus 
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ei{t)^Yo + p'\t),Yo] 

Qp]m= [' A(s,ei(s))ds 

Jo 

Jo I Jo 

and so on. Observe that this procedure preserves the isospectrahty of the flow- 
since the approximation n^'^^{t) lies in so{n) for all m > 1 and t > 0. It is also 
equally possible to develop a formalism based on rooted trees in this case, in 
a similar way as for the standard Magnus expansion. 

Example. The double bracket equation 

Y' = [[F, iV], F], F(0) = Fo e Sym(n) (149) 

was introduced by Brockett [31] and Chu & Driessel [52] to solve certain 
standard problems in applied mathematics, although similar equations also 
appear in the formulation of physical theories such as micromagnctics [179]. 
Here is a constant nxn symmetric matrix. It clearly constitutes an example 
of an isospectral flow with A{t,Y) = [Y,N]. When the procedure (148) is 
applied to (149), one reproduces exactly the expansion obtained in [114] with 
the convergence domain estabhshed in [42]. □ 

3.4 Treatment of general nonlinear equations 

As a matter of fact, the Magnus expansion can be formally generalized to 
any nonlinear explicitly time-dependent differential equation. Given the im- 
portance of the expansion, it has been indeed (re)derived a number of times 
along the years in different settings. We have to mention in this respect the 
work of Agrachev and Gamkrclidzc [3,4,90], and Strichartz [219]. In the con- 
text of Hamiltonian dynamical systems, the expansion was first proposed in 
[191] and subsequently applied in a more general context in [26] with the aim 
of designing new numerical integration algorithms. 

By introducing nonstationary vector fields and fiows, it turns out that one gets 
a linear differential equation in terms of operators which can be analyzed in 
exactly the same way as in section 2. Thus it is in principle possible to build ap- 
proximate solutions of the differential equation which preserve some geometric 
properties of the exact solution. The corresponding Magnus scries expansion 
allows us to write the formal solution and then different approximations can 
be obtained by truncating the series. Obviously, this formal expansion presents 
two difficulties in order to render a useful algorithm in practice: (i) it is not 
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evident what the domain of convergence is, and (ii) some device has to be 
designed to compute the exponential map once the series is truncated. 

Next we briefly summarize the main ideas involved in the procedure. To begin 
with, let us consider the autonomous equation 

x' = f(x), x(0) = xo e M". (150) 

If denotes the exact flow of (150), i.e. x(t) = v9*(xo), then for each inflnitely 
differentiable map g : — > M, g{f^{y)) admits the representation 

g{^\y)) = n9]{y) (151) 

where the operator $* acts on differentiable functions [190]. To be more spe- 
cific, let us introduce the Lie derivative (or Lie operator) associated with f , 

It acts on differentiable functions F : — > M™ (see [6, Chap. 8] for more 
details) as 

LfF(y)=F'(y)f(y), 
where -F'(y) denotes the Jacobian matrix of F. It follows from the chain rule 
that, for the solution </?*(xo) of (150), 

|f((^*(xo)) = (LfF)((^*(xo)), (153) 
and applying the operator iteratively one gets 

^F(^*(xo)) = (Lf^F)(^*(xo)), k>l. 
Therefore, the Taylor series of F((/?*(xo)) at t = is given by 

F((/p*(xo)) = E S(^f^)(^o) = exp(tLf)[F](xo). (154) 

k>o 

Now, putting F{y) — Id(y) = y, the identity map, this is just the Taylor 
series of the solution itself 

<^*(^o) = E T:|(^fId)(xo) = exp(tLf)[Id](xo). 

If we substitute F hj g in (154) and compare with (151), then it is clear that 
$*[5f](y) = exp(tLf)[5(](y). The object exp(tLf) is called the Lie transform 
associated with f . 
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At this point, let us suppose that f (x) can be spht as f (x) = fi(x) + f2(x), in 
such a way that the systems 

X'==fi(x), X' = f2(x) 

have flows and (^2) respectively, so that 

^((^*(y)) = exp(iLfJ[^](y) i = l,2. 

Then, for their composition one has 

gis^l ° V\{y)) = exp(sLfJ exp(iLfJ[^](y). (155) 

This is precisely formula (154) with f = fi, t replaced with s and with F(y) ~ 
exp(tLf2)[(yf](y). Notice that the indices 1 and 2 as well as s and t to the left 
and right of eq. (155) are permuted. In other words, the Lie transforms appear 
in the reverse order to their corresponding maps [99] . 

The Lie derivative Lf satisfies some remarkable properties. Given two functions 
■01, ■02, it can be easily verified that 



Lf (q;i'0i + a2ilJ2) = ai^f^i + a2Lftp2, cti, 0:2 e 

Lf{lpl 1P2) = (Lflpl) 1p2 + Ipl -C'f'02 

and by induction we can prove the Leibniz rule 

i=0 \ i I 



with L\^jJ — Lf (^L\ and Lfi/: = justifying the name of Lie derivative. 
In addition, given two vector fields f and g, then 



aiLf + a2Lg — L^^f+ajg, 

[Lf , Lg] — L{ Lg — Lg Lf = Lh, 

where h is another vector field corresponding to the Lie bracket of the vector 
fields f and g, denoted by h = (f, g), whose components are 

hi^{i,g)i^ LfQi - Lg/, = E (jj ^ - 9 J ^) • ( 156) 

Moreover, from (151) and (153) (replacing F with g) we can write 

i^gK^o) = igiv'i^o)) = (i^f5)(^*(xo)) = $*Lf[^](xo). 
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Particularizing to the function g{-x.) — Idj(x) = Xj, we get 

^$*[Id,](y) = *%(y)[Id,](y), j^l,...,n, y = Xo 

or, in short, 

|$* = $%(y), y = xo, (157) 

i.e., a Unear differential equation for the operator $*. Notice that, as expected, 
equation (157) admits as formal solution 

= exp(iLf(y)), y = Xo. (158) 

We can follow the same steps for the non-autonomous equation 

x' = f(t,x), (159) 

where now the operational equation to be solved is 

= y = Xo. (160) 

To simplify notation, from now on we consider xq as a set of coordinates such 
that f(t,xo) is a differentiable function of xq. Since Lf is a linear operator 
we can then use directly the Magnus series expansion to obtain the formal 
solution of (160) as = exp(Lw(t,xo)), with w = Z^iWj. The first two terms 
are now 



Wi(i,Xo) = / f(s,Xo)ds 
Jo 

W2(t,Xo) = -- / dsi (iS2(f(Sl,Xo),f(s2,Xo)). (161) 

2 Jo Jo 

Observe that the sign of W2 is changed when compared with in (43) and 
the integrals affect only the explicit time-dependent part of the vector field. 
In general, due to the structure of equations (30) and (160), the expression 
for w„(t,xo) can be obtained from the corresponding 0„(t) in the linear case 
by applying the following rules: 

(1) replace A{ti) by f(ti,xo); 

(2) replace the commutator [■, ■] by the Lie bracket (156); 

(3) change the sign in w^(t,xo) for even n. 

Once w'"] = S"=iWi(i,Xo) is computed, there still remains to evaluate the 
action of the Lie transform exp(Lw(t,X(,)) on the initial conditions xq. At time 
t = T, this can be seen as the 1-flow solution of the autonomous differential 
equation 

y' = wN(T,y), y(0)=xo, (162) 
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since y(l) = exp(L^(r,xo))xo = x(T). 

Although this is arguably the most direct way to construct a Magnus expansion 
for arbitrary time dependent nonlinear differential equations, it is by no means 
the only one. In particular, Agrachev and Gamkrelidze [4,90] obtain a similar 
expansion by transforming (160) into the integral equation 

= Id+ f ^'Xsds (163) 
Jo 

which is subsequently solved by successive approximations. Here, for clarity, 
we have denoted Xg = -^f(s,xo)- Then one gets the formal series 

= Id + /* dhXt^ + t dti /*' dt2Xt,Xt, + ■■■ (164) 

JO JO JO 

°° ft ftl ftm-1 

^Id+Y, dh dt2--- dtmXt^---Xt,. 

^1 Jo Jo Jo 

An object with this shape is called a formal chronological series [4], and the 
set of all formal chronological series can be endowed with a real associative 
algebra structure. It is then possible to show that there exists an absolutely 
continuous formal chronological series 

Vt{Xt) dtj dt2--- dtm Gm(Xt, ,...,XtJ (165) 

^iJo Jo Jo 

such that 

Here Gm{Xti, ■ ■ ■ ■, ^t^) are Lie polynomials homogeneous of the first grade in 
each variable, which can be algorithmically constructed. In particular, 

G2{Xtj^,Xt2) 
G3{Xtj^,Xt2,Xt^) 

— * 

The series (165) in general diverges, even if the Lie operator Xt is analytic [3]. 

— * 

Nevertheless, in certain cases convergence holds. For instance, if Xt belongs 
to a Banach Lie algebra B for all t G M, where one has a norm satisfying 
< ||X||||F|| for all X,Y e B and /q \\Xs\\ds = Jq \\L^s,^^^)\\ds < 0.44, 
then Vt{Xt) converges absolutely in B [4]. As a matter of fact, an argument 
analogous to that used in [21] and [172] may allow us to improve this bound 



= Xt, 

= C ([-^t3> [Xt2, XtJ\ + [[Xt3, XjJ, XjJ) 
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and get convergence for 



3.4-1 Treatment of Hamiltonian systems 

We have seen how the algebraic setting we have developed for linear systems 
of differential equations may be extended formally to nonlinear systems. We 
will review next how it can be adapted to the important class of Hamiltonian 
systems. In this context, the role of a Lie bracket of vector fields (156) is played 
by the classical Poisson bracket [226]. 

The Lie algebraic presentation of Hamiltonian systems in Classical Mechanics 
has been approached in different ways and the Magnus expansion invoked in 
this context by diverse authors [162,216,224]. More explicit use of the Magnus 
expansion is done in [191] where the evolution operator for a classical system 
is constructed and its differential equation analyzed. 

To particularize to this situation the preceding general treatment, let us con- 
sider a system with I degrees of freedom and phase space variables x = (q, p) = 

(qi, . . . ,qi,pi, . . . ,pi), where {qi-Pi), i = 1, . . . ,1 are the usual pairs of canoni- 
cal conjugate coordinate and momentum, respectively. By defining the Poisson 
bracket of two scalar functions F(q, p) and G(q, p) of phase space variables 
in the conventional way [226] 



i=l \ 



dFdG dFdG' 
dqi dpi dpi dqi ^ 



we have 



and in particular 



{F,G}- ^ -Q^JiJg^, 



{Xj, Xj} — Jij- 

Here J is the basic symplectic matrix appearing in equation (20) (with n — l)- 
With these definitions the set of (sufficiently smooth) functions on phase space 
acquires the structure of a Lie algebra and we can associate with any such 
function F(x) a Lie operator 

21 Q 

= E g^J.,^. (166) 

which acts on the same set of functions as LpG = {F, G}. It is then a simple 
exercise to show that the set of all Lie operators is also a Lie algebra under 
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the usual commutator [Lp, Lq] — LpLa — LqLf and furthermore 

[Lf: Lc] = L{F,G}- 



Given the Hamiltonian function i7(q, p, t) : M^' x R — >• R, where q, p G R', 
the equations of motion are 

q' = V^H, p' = -V^H, (167) 

or, equivalently, in terms of x, 

x' = JV^H. 

It is then elementary to show that the Lie operator L_h is nothing but the 
Lie derivative Lf (152) associated with the function 

f = JVxi/. 

Therefore, the operational equation (160) becomes 

l$*^ = $*^L_H(y,t), y = xo (168) 

and the previous treatment also holds in this setting. As a result, the Magnus 
expansion reads 

^H^eMLw), (169) 
where W — Zli^i and the first two terms are 



Wii^o)^- tH{xo,s)ds (170) 
Jo 

W^2(xo) = --/ dsi dS2{H{xo,Si),H{xo,S2)}. 

Z Jo Jo 



3.5 Magnus expansion and the Chen-Fliess series 



Suppose that f in equation (159) has the form f(t, x) = J2iLi Ui{t)fi{x) , i.e., 
we are dealing with the nonlinear differential equation 

m 

X'(t) = E^iit) f^(x(^)), X(0) = p, (171) 

i=l 

where Ui{t) are integrable functions of time. Systems of the form (171) appear 
for instance in nonlinear control theory. In that context the functions Ui are 
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the controls and fj arc related to the nonvarying geometry of the system. 
Observe that this problem constitutes the natural (nonlinear) generalization 
of the case studied in section 2.11. 

One of the most basic procedures for obtaining x(T) for a given T is by 
applying simple Picard iteration. For an analytic output function g : — > M, 
from (153) it is clear that 

d 

-gi^m = (L(^„,f^)5)(x(t)) = Y.u,{t){E,g){^{t)), 5(x(0)) = g{p), 

(172) 

where, for simplicity, we have denoted by the Lie derivative Lf.. This can 
be particularized to the case g — Xi, the ith component function. 

By rewriting (172) as an equivalent integral equation and iterating we get 



g{^{t))=g{p)+ / J2^^MiE^,9){Atl))dtl 
= ^(P)+ / E«n(^l) (^«15)(P) + 

/ E««2(^2)(^i2^ny)(x(t2))c^t2Uti (173) 

„t rn I fti ™- / 

+ E uM{Ei^Ei^Ei^g){^{t^))dt^dt2\ dti 



and so on. Notice that in this expression the time dependence of the solution 
is separated from the nonvarying geometry of the system, which is contained 
in the vector fields Ei and need to be computed only once at the beginning 
of the calculation. Next we reverse the names of the integration variables and 
indices used (e.g., rename ii to become is and vice versa), so that 



^(x(i))=^(p)+ E [l^uMdt,^ {Ei,g){p) 

■m / ft Pt2 \ 

+ E E / / uMuMdtM {E,,Ei,g){p) (174) 

m m rn , „t rts rt2 \ 
E E E / / / Ui^{t3)Ui^{t2)Ui^{ti)dtidt2dt3 



^ ^ ^ / ft r^s /■*2 

+ 

{E,,E,,E,,g){p) + 
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Observe that the indices in the Lie derivatives and in the integrals are in 
the opposite order. This procedure can be further iterated, thus yielding the 
formal infinite series 

^(x(i))=^(x(0)) + $^ f ... f' f\,Xts)---Ui,{h)dh-.-dts 

7^in-iJ^ Jo Jo Jo 
Ei^...Ei^g{^{Q)), (175) 

where each ij & L = {1, . . . ,m}. An expression of the form (175) is referred 
to as the Chen-Fliess series, and it can be proved that, under certain circum- 
stances, it actually converges uniformly to the solution of (172) [132]. This 
series originates in K.T. Chen's work [51] on geometric invariants and iterated 
integrals of paths in R". Later, Fhess [87] applied the theory to the analysis 
of control systems. 

One of the great advantages of the Chen-Fliess series is that it can be manipu- 
lated with purely algebraic and combinatorial tools, instead of working directly 
with nested integrals. To emphasize this aspect, observe that each term in the 
series can be identified by a sequence of indices or word w — iii2 • • • in the 
alphabet L through the following two maps: 



M2 ■ w = iii2- ■ -is^ — >(u^l UiXts) i ■■■[ Ui^{ti)dti. . .dts-idt, 

\ Jo Jo Jo 

In fact, the nested integral appearing in the map M.2 can be expressed in a 
simple way, as we did for the linear case in (114) 

ai^...is^ Uisits) •••/ Uij^{ti)dti ■ ■ ■ dts-idts (176) 
Jo Jo Jo 

With this notation, the series of linear differential operators appearing at the 
right-hand side of (175) can be written in the compact form [186] 

E (177) 

where L* denotes the set of words on the alphabet L — {1,2, ...,m}, the 
function a^, is given by (176) for each word w & L* and 

E^^ Ei^-.-Ei^, if w ^ ii- ■ - is e L*. 

It was proved by Chen that the series (177) is an exponential Lie series [51], 
i.e., it can be rewritten as the exponential of a series of vector fields obtained 
as nested commutators of Ei, . . . , E^. Such an expression is referred to in 
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nonlinear control as the formal analogue of a continuous Baker-Campbell- 
Hausdorff formula and also as the logarithm of the Chen-Fliess series [131]. 

Notice the similarities of this procedure with the more general treatment car- 
ried out in subsection 3.4 for the nonlinear differential equation (150). Thus, 
expression (164) constitutes the generalization of (175) to an arbitrary func- 
tion f in (150). Conversely, the logarithm of the Chen-Fliess series can be 
viewed as the corresponding nonlinear Magnus series for the particular non- 
hnear system (171). 

From these considerations, it is clear that, in principle, one can obtain an 
explicit formula for the terms of the logarithm of the Chen-Fliess series in a 
basis of the Lie algebra generated by Ei, . . . , Em, but this problem has been 
only recently solved for any number of operators m and arbitrary order, using 
labelled rooted trees [186]. Thus, for instance, when m = 2, it holds that 

^ a^E^ = eMPiEi + (^^E^ + (5r2[Ei, E2\ + /3ii2[£^i, [Ei, E2\\ 
wei* 

+/?212[^2, [^1,^2]] + •••), (178) 

where, not surprisingly, the expressions of the P coefficients are given by (113) 
with the corresponding change of sign in P12 due to the nonlinear character of 
equation (171). 

Another relevant consequence of the connection between Magnus series and 
the Chen-Fliess series is the following: the Lie series defining the logarithm 
of the Chen-Fliess series can be obtained explicitly from the recurrence (49)- 
(50), valid in principle for the linear case. Of course, the successive terms of 
the Chen-Fliess series itself can be generated by expanding the exponential. 



4 Illustrative examples 

After having reviewed in the preceding two sections the main theoretical as- 
pects of the Magnus expansion and other exponential methods, in this section 
we gather some examples of their application. All of them are standard prob- 
lems of Quantum Mechanics where the exact solution for the evolution oper- 
ator U{t) is well known. Due to their simplicity, higher order computations 
are possible within a reasonable amount of effort. The comparison between 
approximate and exact analytic results may help the reader to grasp the ad- 
vantages as well as the technical difficulties of the methods we have analyzed. 

The examples considered here are treated in [45,136,193,195], although some 
results are unpublished material, in particular those involving highest or- 
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der computations. In subsection 4.1 we present results concerning the most 
straightforward way of dcaUng with ME, namely computations in the Inter- 
action Picture. In subsection 4.2 an application of ME in the adiabatic basis 
is developed. Subsection 4.3 is devoted to illustrate the exponential infinite- 
product expansions of Per and Wilcox. An example on the application of the 
iterative version of ME by Voslamber is given in subsection 4.4. Eventually, 
subsection 4.5 contains an application of the Ploquet-Magnus formalism. 

4-1 ME in the Interaction Picture 

We illustrate the application of ME in the Interaction Picture (see subsec- 
tion 2.9) by means of two simple time-dependent physical systems frequently 
encountered in the literature, for which exact solutions are available: the time- 
dependent forced harmonic oscillator, and a particle of spin ^ in a constant 
magnetic field. In the first case we &x h— 1 for convenience. 

As we will see, ME in the Interaction Picture is appropriate whenever the 
characteristic time scale of the perturbation is shorter than the proper time 
scale of the system. 

To illustrate and evaluate the quality of the various approximations for the 
time-evolution operator, we compute the transition probabilities among non- 
perturbed eigenstates induced by the small perturbation. 

4.1.1 Linearly forced harmonic oscillator 

The Hamiltonian function describing a linearly driven harmonic oscillator 
reads {h — 1) 

H = Ho + V{t), with i/o = ^c^o(p' + g'), V{t)^V2f{t)q (179) 

and f{t) is real. Here q andp stand for the position and momentum operators 
satisfying [q,p] — i and cuq gives the energy level spacing in absence of the 
perturbation V{t). We introduce the usual operators a± = -^{q^ip), so that 
[a_, a+] = 1. With this notation we have 

//o = ^0 (a+a_ + ^) , f{t){a+ + a-). (180) 

The eigenstates of Hq are denoted by \n), so that HqIu) = ujQ{n+ ^)\n) , where 
n stands for the quantum number. With this notation n = corresponds to 
the ground state. 
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For simplicity in the computations we choose ujq = 1. In accordance with the 
prescriptions in section 2.9, the Hamiltonian in the Interaction Picture is given 
by (92) and reads 

Hrit) = e^-^°* V{t) e-^^°* = f{t){e^* a+ + e"^* a_). (181) 

Accordingly, the evolution operator is factor ized as 

U{t, 0) = exp{-iHot) Ui{t, 0), (182) 

where the new evolution operator Uj is obtained from Uj — Hi{t)Ui = 
-iHi{t)Ui. 

The infinite Magnus series terminates in the present example. It happens be- 
cause the second order Magnus approximant, which involves the computation 
of 

[Hj{h), Hj{h)] = fih)fih) (e^(*^-*^) [a+, a_] + e''^'^-'^^ [a_, a+]) 

= 2i/(ti)/(t2)sin(t2-ti) (183) 

reduces to a scalar function. Thus Magnus series in the Interaction Picture 
furnishes the exact evolution operator irrespective of f{t): 

= exp — i{aa^ + a*a_) — ifi^ (184) 
= exp(— iQ;a+) exp(— iQ;*a_) exp(— — |q;|^/2), 

where we have defined 

a= ['dhf(h)e''\ (185) 
Jo 

P= tdt, rdt2f{h)f{t2)sm{t2-t,). (186) 
Jo Jo 

Equations (184) and (182) yield the exact time-evolution operator for the 
linearly forced harmonic oscillator Hamiltonian (179) [140]. 

To compute transition probabilities between free harmonic oscillator states of 
quantum numbers n and m, 

Pn-.m^\{m\Ui\n)\'', (187) 

the last form in (184) results most convenient. Specifically, assuming that the 
oscillator was initially in its ground state |0), we get in particular the familiar 
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Poisson distribution for the transition probabilities 



Po^n = -r kl^'' exp(-|a|2). 
nl 



(188) 



4- 1-2 Two-level quantum systems 

The generic Hamiltonian for a two-level quantum system can be written down 
in the form 



H{t)^ 



^ E^{t) Cit)^ 



(189) 



\C*it) E,it)J 

where Ei{t), E2{t) are real functions and C{t) is, in general, a complex function 
of t. We define the solvable piece of the Hamiltonian as the diagonal matrix 



Ho{t) 



Ei{t) 
E^it) 



(190) 



and all the time- dependent interaction described by the function C{t) is con- 
sidered as a perturbation. In the Interaction Picture the new Hamiltonian 
reads (see (92)) 



/ 



Hi{t) 







C*(t)exp (^-i dt'uit') 



C{t)exp(^i J^dt'u{t') 







(191) 



with uj = {El —E2)/h. Suppose now that Hq is time-independent. Then U{t) ~ 
exp{HQt) Uj{t). Without loss of generality, the Hq may be rendered traceless, 
so that El = —E2 = E. Thus ±E denote the eigenenergies associated to the 
eigenvectors |-|-) = (1,0)^, |— ) = (0,1)^ of Hq, the unperturbed system. In 
terms of Pauli matrices the Hamiltonian in this case may be expressed as 



H{t)^-nu;as + f{t)ai + g{t)a2, 



(192) 



where / = Re(C) and g = — Im(C). 



Since Hq is diagonal, the transition probability between eigenstates |-|-), |— ) 
of Hq is simply 

P{t) = \{+\Uj{t)\-)f. (193) 
As the evaluation of (193) requires the computation and manipulation of ex- 
ponential matrices involving Pauli matrices, formulas (18) and (19) in section 
1.2 come in hand here. 

Next, we study two particular cases of interaction for which the exact solution 
of the time evolution operator admits an analytic expression. 
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1- Rectangulair step. Suppose that in (192) g — 0, namely, 

H{t) = ^hu;a3 + f{t)ai (194) 

with / = for t < and / = Vq for t > 0. Ahernatively, if we restrict ourselves 
to compute an observable such as the transition probability, this example is 
equivalent to a rectangular mound (or rectangular barrier) of width T — t 
above. The exact solution for this problem reads 

U{t, 0) = exp (-Z (|a3 + ^ai) , (195) 
which yields the exact transition probability 

^e.= 4:^sinV7^ + e/4 (196) 

between eigenstates |+), |— ) of Hq. Here we have denoted 7 = Vot/h and 
^ = u;t. 

The Interaction Picture is defined here by the exphcit integration of the diag- 
onal piece in the Hamiltonian, so that U — exp{—i^a3/2)Ur, where Uj stands 
for the time evolution operator in the Interaction Picture and obeys 

U'j = Hi{t) Ui, Uj{0) = I (197) 

with 

Hi{t) = f{t){aiCos^-a2sm^). (198) 

A computation with the usual time-dependent perturbation theory gives for 
the first orders (formula (52)) 



4'^=4'^ = ^sin2(e/2) (199) 



p(3) ^ p(4) ^ 7^ 
Pt Pt C2 




9 sin - -I- sin ^ — 6f cos - -|- 4 sin^ - 
2 2^2 2, 



2 



Notice that Ppt > 1 may happen in the equations above because the unitary 
character of the operator U (t) is not preserved by the usual time-dependent 
perturbation formalism. 

In this example it is not difficult to compute the first four terms in the Magnus 
series corresponding to Uj{t) = expQ(t) in (197). To facilitate the notation, 
we define s = sin^ and c = cos^. The Magnus approximants in the Interaction 
Picture may be written down in terms of Pauli matrices and read explicitly 
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^4 



^^3 



^2 



I 



I 




3 



(c + 8)s]. 



{5 + c)s]+a2m-s)s-^l-c)]} 



(200) 



The first two formulae for the approximate transition probabihties are, respec- 
tively 



PS^--k-^sm\C/2), A=[4sin^(e/2) + ^(sine-0T/^- (201) 



We omit explicit expressions for P^-^ and P^j' since they are quite involved. 
However, we include their outputs in Figures 3, 4 and 5, where we plot the first 
to fourth order approximate transition probabilities with ME in the Interaction 
Picture and compare them to the exact case and also with perturbation theory 
outputs. In Figure 3 and 4 we set 7 = 1.5 and 7 = 2 respectively, whereas in 
Figure 5 we fix ^ = 1. 

Wc observe that for the Magnus expansion in the Interaction Picture, the 
smaller the value of the parameter ^ the better works the approximate solution. 
As a matter of fact, in the sudden hmit, ^ <S 1, ME furnishes the exact 
result; unlike perturbation theory. As far as the intensity of the perturbation 
7 increases, the quality of the approximations spoils. This effect is much more 
dramatic for the standard perturbation theory. 

On the other hand, it is clear from (198) that 



whence \\Hj(ti)\\2 dti = 7, and thus Theorem 9 guarantees that the Mag- 
nus expansion in the Interaction Picture is convergent if 7 < tt. Notice that 
this is always the case for the parameters considered in Figures 3-5. The esti- 
mate 7 < TT for the convergence domain in the Interaction Picture should be 
compared with the corresponding one in the Schrodinger picture: 




—00 



Hj{h)\\2dt,^ I \f{h)\dh^Vot, 
Jo 
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Fig. 3. Rectangular step: Transition probabilities as a function of ^, with 7 = 1.5. 
The solid line corresponds to the exact result (196). Broken lines stand for ap- 
proximations obtained via ME and lines with symbols correspond to perturbation 
theory, according to the legend. Computations up to fourth order, in the Interaction 
Picture. 

Notice then that, as pointed out in subsection 2.9, a change of picture allows 
us to improve the convergence of the Magnus expansion. 

2- Hyperbolic secant step: Rosen— Zener model. In the Rosen Zener 
Hamiltonian [204] the interaction C{t) in (189) is given by the real function 
V{t) = Vosech(t/T), where T determines the time-scale. We will use the 
notation 7 = irVoT/h and ^ = uT = 2ET/h. 

The corresponding Hamiltonian in terms of Pauli matrices is 

H{t) = Eas + V{t)ai = a{t) ■ a, V{t) = Vq/ cosh(t/T), (202) 

with a = {V{t), 0, E). In the Interaction Picture one has 

Hj{s) = V{s){ai cos(^s) - a2 sin(^s)) (203) 

in terms of the dimensionless time- variable s = t/T. Notice that ^ measures 
the ratio between the interaction time T and the internal time of the system 
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Fig. 4. Rectangular step: Transition probabilities as a function of ^, with 7 = 2. 
Lines are coded as in Figure 3. Computations up to fourth order, in the Interaction 
Picture. 

h/2E. From (203), and after straightforward calculation, the first and second 
ME operators are readily obtained. 

The exact result for the transition probability (provided the time interval ex- 
tends from —00 to +00), as well as perturbation theory and Magnus expansion 
up to second order read [195] 
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Fig. 5. Rectangular step: Transition probabilities as a function of 7, with ^ = 1. 
Lines are coded as in Figure 3. Computations up to fourth order, in the Interaction 
Picture. 



In Figures 6 and 7 we plot some results from the formulae in (204). In Figure 
6 we take 7 = 1.5 and in Figure 7 we set ^ = 0.3. 

Similarly to the case of the rectangular step, we observe in Figure 6 that the 
Magnus expansion works better in the sudden regime defined by ^ <^ 1, namely, 
when the internal time of the system h/2E is much larger than the time scale 
T of the perturbation. Also, Figure 7 illustrates how the approximations spoil 
as far as the intensity 7 increases. Notice in both figures the unitarity violation 
of the approximation built with the usual perturbation theory. 

In this case a simple calculation shows that 

\\Hi{t)hdt = {l/h) / \V{t)\dt = Vo7rT/h = j, 

-00 J 00 

and thus the Magnus series converges at least for 7 < tt. 
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Fig. 6. Rosen-Zener model: Transition probabilities (204) as a function of S^, with 
7 = 1.5. The solid line stand for the exact result. Broken lines stand for approxima- 
tions obtained via ME and triangles correspond to perturbation theory, according 
to the legend. Computations up to second order, in the Interaction Picture. 

4-2 ME in the Adiabatic Picture 



Here we will illustrate the effect of using the Adiabatic Picture introduced in 
Section 2.9. The use of this type of preliminary transformation is convenient 
whenever the time scale of the interaction is much larger than the proper time 
of the unperturbed system. 

Since an adiabatic regime conveys a smooth profile for the perturbation, 
namely, existence of derivatives, the case of the rectangular step cannot be 
properly used for the sake of illustration. 



4-2.1 Linearly forced harmonic oscillator 

For the linearly driven harmonic oscillator the procedure yields the exact so- 
lution, as in the preceding subsection 4.1.2, albeit the method is a bit more 
involved technically [140]. 
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Fig. 7. Rosen-Zener model: Transition probabilities as a function of 7, with = 0.3. 
Lines are coded as in Figure 6. Computations up to second order, in the Interaction 
Picture. 

4-2.2 Rosen-Zener model 

Following [138] we will deal with the Rosen-Zener model (see Section 4.1.2 
and [195]) since it allows a clear illustration of the adiabatic regime. 

The preliminary linear transformation G{s) defined in (94) for the Hamiltonian 
(202) is given by G{s) = b ■ cr, where the unit vector b points in the direction 
b = a + k [k =unit vector along the z-axis). Remind that s = t/T is the 
dimensionless time-variable. The evolution operator gets then factorized as 

Ug{s)=G\s)U{s)G{so), (205) 
which according to (90), satisfies the equation 

dUc 



ds 

Here Hq = —iHc/h is given by 



Hg{s)Ug. (206) 



Hg{s) = ^aa^ - t^a2, (207) 
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with — Eq + Vq/ cosh^ s and cot 9 — (Eo/Vo) cosh s. 

Next, in analogy to (203), we introduce the Adiabatic Interaction Picture 
which allows us to integrate the diagonal piece of Hg{s)- The time-evolution 
operator gets eventually factorized as 

Ug{s) = exp [{-iT/h) ds'a{s')a^ U^\s)exp [{-iT/h) ds'a(s>3^ , 

(208) 



where Uq\s) obeys the equation 



^^H^J\s)uil\ (209) 



with 
and 



ds 

H^J^ (s) = -i(e'/2) [ai sin A(s) + (72 cos A{s)\ , (210) 

Ais) = ^ /' ds' ais') = I In \^ + ^ arctan ^p. (211) 
Ti Jo 2 1 — p TT TTt 



We have introduced the definition 

p = {l-[l + (7r{/27)]sin2^}V2, (212) 

in terms of the dimensionless strength parameter 7 = nVoT/h and 6. Using 
the ME to first order in the adiabatic basis (which coincides with the fixed 
one at s = ±00) one finds the spin-fiip approximate (first order) transition 
probability 

dOsm A{s{e)) 



'^'^Pj]) = sm^ 



(213) 



In Figure 8 we compare the numerical results given by the new approximation 
(213) with the exact formula Pgx in (204). For the sake of illustration we plot 
also the results in the usual Interaction Picture to first order in ME (see Pj^^ 
in (204)). It is noteworthy the gain achieved when using the adiabatic ME in 
the intermediate regime, namely, moderate values of ^, although only the first 
order is considered. Here the adiabatic regime corresponds to large values of 
{ = 2ET/h {e = 1/T < 1). 



It should be also noticed that for the Hamiltonian (210) one has 
£ \\H^J\s,)hds, = ^|^(.) - ^(.o)| < ^27r = TT 



and thus the convergence condition given by Theorem 9 is always satisfied. In 
other words, for this example the Magnus expansion is always convergent in 
the Adiabatic Interaction Picture. 

More involved illustrative examples of ME in the Adiabatic Picture may be 
found in [139,140]. 
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Fig. 8. Rosen-Zener model: Transition probability as a function of ^, with 7 = 1. The 
solid line stands for the exact result in (204). The remaining lines stand for first 
order computations with ME (circles) and perturbation theory (triangles). Open 
symbols are for the Adiabatic Picture and solid symbols for the Interaction Picture. 

4-3 Fer and Wilcox infinite-product expansions 



Next we illustrate the use of Fer and Wilcox infinite-product expansions by 
using the same time- dependent systems as before. 



4-3.1 Linearly forced harmonic oscillator 

Since the commutator [a_ , a+j is a c-number, Fer iterated Hamiltonians H^"^^ 
with n > 1 eventually vanish so that F„ = for n > 2. The Wilcox operators 
Wn with n > 2 in eq.(106) vanish for the same reason. Thus, in this particular 
case, the second-order approximation in either method leads to the exact so- 
lution of the Schrodinger equation just as ME did. To sum up, the final result 
reads 

JJ^ ^ gHi+Ha ^ gWi ^W2 ^ gFi gF2 ^ g-i(aa++a*a_) ^ ^214) 

where a(t) and P(t) are given in (185) and (186), respectively. 
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4-3.2 Two-level quantum system: Rectangular step 

For the Hamiltonian (194), the first-order Fer and Wilcox operators in the 
Interaction Picture verify 

F.^Wi^Qi^ f'dtiHiiti), (215) 
Jo 

where Hj{t) is given by (198). The exphcit expression is collected in (200) 
(first equation). Analogously, the second equation there also corresponds to 
the second-order Wilcox operator W2 — 0,2- 

To proceed further with Fer's method we must calculate the modified Hamil- 
tonian H^^^ in (96). After straightforward algebra one eventually obtains 

where 9 = (27/^) sin(,^/2) (notice that H^^^ and therefore F2 depend on cti and 
(72, while W2 is proportional to as). Since it does not seem possible to derive 
an analytical expression for F2, the corresponding matrix elements have been 
computed by replacing the integral by a conveniently chosen quadrature. 

The transition probability P{t) from an initial state with spin up to a state 
with spin down (or viceversa) is given by (193). This expression has been 
computed on assuming: Uj ^ e^i = e^^, Ui ~ e^^e^'^^ Uj ~ c^^^ jj^ ^ 
e^i e^2 e-^-^ and [// ^ e^'-^ e^-^, and the results have been compared with the 
exact analytical solution (196). 

In Figures 9 and 10 we show the transition probability P as a function of ^ 
for two different values of 7 (7 = 1.5 and 7 = 2, respectively), while in Figure 
11 we have plotted P versus 7 for ^ fixed. Notice that the second order in 
the Wilcox expansion does not contribute to the transition probability (this 
is similar to what happens in perturbation theory). On the other hand, Fer's 
second-order approximation is already in remarkable agreement with the exact 
result whereas the third order cannot even be distinguished from the exact 
result in Figure 10 at the cost of a much larger computational effort. Wilcox 
approximants preserve unitarity but do not give acceptable approximations. 

4.4 Voslamber iterative method 

Just to keep the same structure as in preceding subsections, we mention that 
the Voslamber iterative method of subsection 3.1 also yields the exact solu- 
tion for the linearly driven harmonic oscillator after computing the second 
iteration. 
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Fig. 9. Rectangular step: Transition probability as a function of ^, with 7 = 1.5. 
The solid line corresponds to the exact result (196). Broken lines stand for Fer and 
Wilcox approximations up to third order. Computations up to third order, in the 
Interaction Picture. 

Next, for the two- level system with a rectangular step described by the Hamil- 
tonian (194) we compute the second iterate and compare with second or- 
der ME approximation for Uj. As a test, we shall obtain again the transition 
probability P{t) given by (196). The expression (193) will be calculated here 
on assuming: Uj ~ exp il^^^ = exp [// ~ exp(f2i + ^2) and Uj ~ exp fi*^^-*. 

The second order Magnus approximation to the transition probability is given 
by (201), whereas the second iterate is obtained from (124), 

f](2)(t) = -il r*{[sin2(A) + cos2(A) cose] ^1 

U! Jo 

+ cos^(A)sinea2 -sin(2A)sin(e/2)(T3}dC, (217) 

where A = ^|sin(^/2)|, 7 = Vot/h, ^ = ut. Since it does not seem possible 
to derive an analytical expression for fi'-^-' , the corresponding matrix elements 
have been computed by approximating the integral in (217) with a sufficiently 
accurate quadrature. 
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Fig. 10. Rectangular step: Transition probability as a function of with 7 = 2. 
Lines are coded as in Figure 9. Computations up to third order, in the Interaction 
Picture. 



In Figure 12, the various approximated transition probabilities as well as the 
exact result (196) have been plotted as a function of ^ for a fixed value 7 = 1.5. 
We observe that the approximation from the second iterate keeps the trend of 
the exact solution in a better way that the second order Magnus approximation 
does. 

In Figure 13 we have plotted the corresponding transition probabilities versus 
7 for fixed value ^ = 1. Although locally the second order Magnus approxi- 
mation may be more accurate, it seems that the trend of the exact solution is 
mimicked for a longer interval of 7. 

As it has already been pointed out above, the Magnus expansion works the 
better the more sudden the perturbation. Thus, the re-summation involved in 
the iterative method improves a bit that issue. Further results on the present 
example may be found in [193]. 
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Fig. 11. Rectangular step: Transition probability as a function of 7, with ^ = 1. 
Lines are coded as in Figure 9. Computations up to third order, in the Interaction 
Picture. 

4-5 Linear periodic systems: Floquet-Magnus formalism 



Next we deal with a periodically driven harmonic oscillator, where the infi- 
nite expansions obtained from Floquet-Magnus recurrences (see subsection 
3.2), are utterly summed. Comparison with the exact solution illustrates the 
feasibility of the method. 

The particular system we consider is described by the Hamiltonian (179) with 
fit) = ^ cos cut and ujq < 00. Once the recurrences of the Floquet-Magnus 
expansion (see section 3.2) are explicitly computed for several orders, their 
general term may be guessed by inspection. For the Floquet operator we get 




and the associated transformation results from 
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Fig. 12. Rectangular step: Transition probability in the two level system as a func- 
tion of ^, for 7 = 1.3: Exact result (204) (solid line), second iterate of Voslamber 
method (dashed line), second order ME (dotted line) and first Voslamber iterate (or 
order in ME) (dash-dotted line). Computations are done in the Interaction Picture. 



s,„M)£(. + i)f^)'-ff (^)' 



fc=0 



fe=0 



i [sin (out) q + ujq (cos (out) — 1) p] — X! ( — 



(219) 



The resulting series may be summed in closed form thus yielding the Floquet 
operator 



F = -i 



P 



— z 



Auo (p2 - 1) 



(220) 



with p = uj/ujq. Its eigenvalues are the so-called Floquet eigenenergies [144] 
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Fig. 13. Rectangular step: Transition probability in the two level system as a func- 
tion of 7, for ^ = 1. Lines are coded as in Figure 12. Computations are done in the 
Interaction Picture. 

The corresponding A transformation after summation of the series in (219) is 



1 -p2 



(p sin ujt) q + (cos ut — l)p + 



2p2 



1-p^ 



+ cos ut 



P sin ut 



(222) 

Notice that, as they should, both operators are skew-Hermitian and reproduce 
the exact solution of the problem. 



5 Numerical integration methods based on the Magnus expansion 



5. 1 Introduction 



The Magnus expansion as formulated in section 2 has found extensive use in 
mathematical physics, quantum chemistry, control theory, etc, essentially as 
a perturbative tool in the treatment of the linear equation 

Y' = A{t)Y, Y{to) = Yo. (223) 
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When the recurrence (49)- (50) is apphed, one is able to get expUcitly the suc- 
cessive terms flk in the series defining Q as hncar combinations of muhivariate 
integrals containing commutators acting iteratively on the coefficient matrix 
A, as in (42)-(46). As a result, with this scheme analytical approximations 
to the exact solution are constructed explicitly. These approximate solutions 
are fairly accurate inside the convergence domain, especially when high or- 
der terms in the Magnus series are taken into account, as illustrated by the 
examples considered in section 4. 

There are several drawbacks, however, involved in the procedure developed so 
far, especially when one tries to find accurate approximations to the solution 
for very long times. The first one is implicitly contained in the analysis done 
in section 2: the size of the convergence domain of the Magnus series may be 
relatively small. The logarithm of the exact solution Y{t) may have complex 
singularities and this implies that no series expansion can converge beyond 
the first singularity. This disadvantage may be up to some point avoided by 
using different pictures, e.g. the transformations shown in section 2.9, in or- 
der to increase the convergence domain of the Magnus expansion, a fact also 
illustrated by several examples in section 4. Unfortunately, these preliminary 
transformations sometimes either do not guarantee convergence or the rate 
of convergence of the series is very slow. In that case, accurate results can 
only be obtained provided a large number of terms in the series are taken into 
account. 

The second drawback is the increasingly complex structure of the terms Q,k in 
the Magnus series: each Qk is a ^-multivariate integral involving a linear com- 
bination of (k — l)-nestcd commutators of A evaluated at different times tj, 
i = 1, . . . ,k. Although in some cases these expressions can be computed explic- 
itly (for instance, when the elements of A and its commutators are polynomial 
or trigonometric functions), in general a special procedure has to be designed 
to approximate multivariate integrals and reduce the number of commutators 
involved. 

When the entries of the coefficient matrix A{t) are complicated functions of 
time or they are only known for certain values of t, numerical approximation 
schemes are unavoidable. In many cases it is thus desirable to obtain just 
numerical approximations to the exact solution at many different times. This 
section is devoted precisely to the Magnus series expansion as a tool to build 
numerical integrators for equation (223). 

Before embarking ourselves in exposing the technical details contained in this 
construction, let us first introduce several concepts which are commonplace in 
the context of the numerical integration of differential equations. 
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Given the general (nonlinear) ordinary differential equation (ODE) 

x' = f(t,x), x(to)=xoeC^ (224) 



standard numerical integrators, such as Runge-Kutta and multistep methods, 
proceed as follows. First the whole time interval [to,t/], is split into subin- 
tcrvals, [tn^i,tn], n = 1,. . . ,N, with tjy = tf, and subsequently the value of 
x(t„) is approximated with a time-stepping advance procedure of the form 



starting from xq. Here the map $ depends on the specific numerical method 
and hn — tn+i — tn are the time steps. For simplicity in the presentation we 
consider a constant time step h, so that tn — to + nh. In this way one gets 
x„+i as an approximation to x(t„_|_i). In other words, the exact evolution of the 
system (224) is replaced by the discrete or numerical flow (225). The simplest 
of all numerical methods for (224) is the explicit Euler scheme 



It computes approximations x„ to the values x(t„) of the solution using one 
explicit evaluation of f at the already computed value x„_i. In general, the 
numerical method (225) is said to be of order p if, assuming x„ = x(i„), then 
x^+i = x(i„+i) + 0{h^^^). Thus, in particular, Euler method is of order one. 

Of course, more elaborate and efficient general purpose algorithms, using sev- 
eral f evaluations per step, have been proposed along the years for the numer- 
ical treatment of equation (224). In fact, any standard software package and 
program library contains dozens of routines aimed to provide numerical ap- 
proximations with several degrees of accuracy, including (explicit and implicit) 
Runge-Kutta methods, linear multistep methods, extrapolation schemes, etc, 
with fixed or adaptive step size. They are designed in such a way that the 
user has to provide only the initial condition and the function f to obtain 
approximations at any given time. 

This being the case, one could ask the following question: if general purpose 
integrators are widely available for the integration of the linear equation (223) 
(which is a particular case of (224)), what is the point in designing new and 
somehow sophisticated algorithms for this specific problem? 

It turns out that in the same way as for classical time-dependent perturbation 
theory, the qualitative properties of the exact solution are not preserved by the 

numerical approximations obtained by standard integrators. This motivates 
the study of the Magnus expansion with the ultimate goal of constructing 
numerical integration methods. We will show that highly accurate schemes 
can be indeed designed, which in addition preserve qualitative properties of 



x„+i = $(/l„, X, 



. . . ,/io,Xo) 



(225) 




(226) 
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the system. The procedure can also incorporate the tools developed for the 
analytical treatment, such as preliminary linear transformations, to end up 
with improved numerical algorithms. 

We first summarize the main features of the well know class of Runge-Kutta 
methods as representative of integrators of the form (225). They are introduced 
for the general nonlinear ODE (224) and subsequently adapted to the linear 
case (223). 



5.2 Runge-Kutta methods 



The Runge-Kutta (RK) class of methods are possibly the most frequently used 
algorithms for numerically solving ODEs. Among them, perhaps the most suc- 
cessful during more than half a century has been the 4th-order method, which 
applied to equation (224) provides the following numerical approximation for 
the integration step i„ i— > tn+i — tn + h: 



Y2 = x„ + ^f(t,,Yi) 

Y3 = x„ + ^f(i„ + ^,Y2) (227) 

h 

Y4 = X„ + /lf(t„ + -,Y3) 
Xn+l = Xn + ^ ( f (t„, Yi) + 2i{tn + ^, Y2) + 2f (t, + ^, Y3) + f (t, + h, Y4) ] . 



Notice that the function f can always be computed explicitly because each 
Yj depends only on the Yj, j < i, previously evaluated. To measure the 
computational cost of the method, it is usual to consider that the evaluation 
of the function f (t, x) is the most consuming part. In this sense, scheme (227) 
requires four evaluations, which is precisely the number of stages (or inner 
steps) in the algorithm. 

The general class of s-stage Runge-Kutta methods are characterized by the 
real numbers aij, bi {i, j — 1, . . . , s) and Cj = Sj=i <^ij7 ^ 

s 

Yj = x„ + /?. ^ aij f{tn + Cjh, Yj), i = 1, . . . , s 

s 

X-n+l — X^ + hJ2bii{tn + Cih,^i), (228) 

i=l 
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where Yj, i = l,...,s are the intermediate stages. For simphcity, the as- 
sociated coefficients are usuaUy displayed with the so-called Butcher tableau 
[37,100] as follows: 



Cl 


ail ■ 


■ ais 




dsl ■ 






bi . 


■ bs 



(229) 



If aij = 0, j > then the intermediate stages Yj can be evaluated recursively 
and the method is explicit. In that case the zero aij coefficients (in the upper 
triangular part of the tableau) are omitted for clarity. With this notation, 'the' 
4th-order Runge-Kutta method (227) can be expressed as 
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Otherwise, the scheme is implicit and requires to numerically solve a system 
of s d nonlinear equations of the form 



y = X„ + /iG(/i,y), (231) 

where y = (Yi, . . . , Y^)^, X^j = (x„,...,x„)^ e W^, and G is a function 
which depends on the method. A standard procedure to get x„+i from (231) 
is applying simple iteration: 

y[^l - X„ + G(/i, y[^-il), J = 1, 2, . . . (232) 

When h is sufficiently small, the iteration starts with y[°l = X„ and stops once 
||yb]_yb-i] II is smaller than a prefixed tolerance. Of course, more sophisticated 
techniques can be used [100] . 

After these general considerations, let us turn our attention now to the linear 
equation (223). When dealing with numerical methods applied to this equa- 
tion, it is important to keep in mind that the relevant small parameter here is 
no longer the norm of the matrix A{t) as in the analytical treatment, but the 
time step h. For this reason, the concept of order of accuracy of a method is 
different in this context. Now f(t,x) = /l(t)x and the general class of s-stage 
Runge-Kutta methods adopt the more compact form 
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i=l 

with = A{tn + Cj/i). In terms of matrices, this is equivalent to 



(233) 
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(234) 



so that the apphcation of the method for the integration step t„ i— > 
tn + h can also be written as 



x„+i = (ld + h( biAi ■ ■ ■ bsAs j (isd - hA^ ^ Isdxdj Xn, 



(235) 



where Igdxd = {Idi ^d, ■ ■ ■ , Id)^ and is the dxd identity matrix. For instance, 
taking s = 2 and using Matlab this can be easily implemented as follows 



A^[all*Al al2 * A2 ] a2l * Al a22*^2]; 

x^{Id + h*[hl*Al h2*A2]*{{nd-h*A)\[Id Id]'))*x; 

where bl, 62, all, . . . , a22 are the coefficients of the method, Al, A2 correspond 
to Al, A2 and Id, I2d are the identity matrices of dimension d and 2d, respec- 
tively. 

There exists an extensive literature about RK methods built for many different 
purposes [37,100]. It is therefore reasonable to look for the most appropriate 
scheme to be used for each problem. In practice, explicit RK methods are pre- 
ferred, because its implementation is usually simpler. They typically require 
more stages than implicit methods and, in general, a higher number of evalu- 
ations of the matrix A{t), although this is not always the case. For instance, 
the 4-stagc fourth order method (227) only requires two new evaluations of 
A(t) per step (A(t„ + h/2) and A(tn + h)), which is precisely the minimum 
number of evaluations needed by any fourth order method. In our numerical 
examples we will also use the 7-stage sixth order method with coefficients [37, 
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Observe that this method only requires three new evaluations of the matrix 
A{t) per step. This is, in fact, the minimum number for a sixth-order method. 
Other implicit RK schemes widely used in the literature involving the mini- 
mum number of stages at each order are based on Gauss-Legendre collocation 
points [100]. For instance, the corresponding methods of order four and six 
(with two and three stages, respectively) have the following coefficients: 
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5.3 Preservation of qualitative properties 

Notice that the numerical solution provided by the class of Runge-Kutta 
schemes, and in general by an integrator of the form (225), is constructed 

as a sum of vectors in W^. Let us point out some (undesirable) consequences 
of this fact. Suppose, for instance, that x is a vector known to evolve on a 
sphere. One does not expect that x„ built as a sum of two vectors as in (225) 
preserves this feature of the exact solution, whereas approximations of the 
form x„+i = Qn^n, with Qn an orthogonal matrix, clearly lie on the sphere. 

In section 3.3 we have introduced isospectral flows (143), which include as a 
particular case the system 

Y'^[A{t),Y], Y{to)^Yo, (238) 

with A a skew-symmetric matrix. As we have shown there, if Yq is a symmetric 
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matrix, the exact solution can be factorized as Y(t) = Q(t)YoQ'^(t), with Q{t) 
an orthogonal matrix satisfying the equation 

Q' = A{t)Q, g(0) = I (239) 

and, in addition, Y{t) and Y{0) have the same eigenvalues. This is also true 
when Y is Hermitian, A is skew-Hermitian and Q is unitary, in which case 
(238) and (239) can be interpreted as particular examples of the Heisenberg 
and Schrodinger equations in Quantum Mechanics, respectively. 

When a numerical scheme of the form (225) is applied to (239), in general, the 
approximations Qn will no longer be unitary matrices and therefore Y^ and 
Yq will not be unitarily similar. As a result, the isospectral character of the 
system (238) is lost in the numerical description. Observe that explicit Runge- 
Kutta methods employ the ansatz that locally the solution of the differential 
equation behaves like a polynomial in t, so that one cannot expect that the 
approximate solution be a unitary matrix. In this sense, explicit Runge-Kutta 
methods present the same drawbacks in the numerical analysis of differential 
equations as the standard time-dependent perturbation theory when looking 
for analytical approximations. 

Implicit Runge-Kutta methods, on the other hand, can be considered as ratio- 
nal approximations, and in some cases the outcome they provide is a unitary 
matrix. For this class of methods, however, matrices of relatively large sizes 
have to be inverted, making the algorithms computationally expensive. Fur- 
thermore, the evolution of many systems (including highly oscillatory prob- 
lems) can not be efficiently approximated neither by polynomial nor rational 
approximations. 

With all these considerations in mind, a pair of questions arise in a quite 
natural way: 

(Ql) Is it possible to design numerical integration methods for equation (224) 
such that the corresponding numerical approximations still preserve the 
main qualitative features of the exact solution? 

(Q2) Since the Magnus expansion constitutes an extremely useful procedure for 
obtaining analytical (as opposed to numerical) approximate solutions to the 
linear evolution equation (223), is it feasible to construct efficient numerical 
integration schemes from the general formulation exposed in section 2? 

It turns out that both questions can be answered in the affirmative. As a 
matter of fact, it has been trying to address (Ql) how the field of geomet- 
ric numerical integration has been developed during the last years. Here the 
aim is to reproduce the qualitative features of the solution of the differen- 
tial equation which is being discretised, in particular its geometric properties. 
The motivation for developing such structure-preserving algorithms arises in- 
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dependently in areas of research as diverse as celestial mechanics, molecular 
dynamics, control theory, particle accelerators physics, and numerical analysis 
[99,119,146,165,166]. 

Although diverse, the systems appearing in these areas have one important 
common feature. They all preserve some underlying geometric structure which 
influences the qualitative nature of the phenomena they produce. In the field 
of geometric numerical integration these properties are built into the numer- 
ical method, which gives the method an improved qualitative behaviour, but 
also allows for a significantly more accurate long-time integration than with 
general-purpose methods. 

In addition to the construction of new numerical algorithms, an important 
aspect of geometric integration is the explanation of the relationship between 
preservation of the geometric properties of a numerical method and the ob- 
served favourable error propagation in long-time integration. 

Geometric numerical integration has been an active and interdisciplinary re- 
search area during the last decade, and nowadays is the subject of intensive 
development. Perhaps the most familiar examples of geometric integrators are 
symplectic integration algorithms in classical Hamiltonian dynamics, symmet- 
ric integrators for reversible systems and methods preserving first integrals and 
numerical methods on manifolds [99]. 

A particular class of geometric numerical integrators are the so-called Lie- 
group methods. If the matrix A{t) in (223) belongs to a Lie algebra q, the 
aim of Lie-group methods is to construct numerical solutions staying in the 
corresponding Lie group Q [119]. 

With respect to question (Q2) above, it will be the subject of the next subsec- 
tion, where the main issues involved in the construction of numerical integra- 
tors based on the Magnus expansion are analyzed. The methods thus obtained 
preserve the qualitative properties of the system and in addition are highly 
competitive with other more conventional numerical schemes with respect to 
accuracy and computational effort. They constitute, therefore, clear examples 
of Lie-group methods. 

5.4 Magnus integrators for linear systems 

Since the Magnus series only converges locally, as we have pointed out before, 
when the length of the time-integration interval exceeds the bound provided 
by Theorem 9, and in the spirit of any numerical integration method, the usual 
procedure consists in dividing the time interval [to,^/] into N steps such that 
the Magnus series converges in each subinterval [tn-i, tn], n — 1, . . . ,N, with 
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tN — tf. In this way the solution of (223) at the final time tf is represented by 

N 

YM = n cxp(n(i„, tn-i)) Yo, (240) 

n=l 

and the series 0(t„,t„_i) has to be appropriately truncated. 

Early steps in this approach were taken in [49] and [50], where second and 
fourth order numerical schemes were used for calculations of coUisional inelas- 
ticity and potential scattering, respectively. Those authors were well aware 
that the resulting integration method "would become practical only when the 
advantage of being able to use bigger step sizes outweighs the disadvantage in 
having to evaluate the integrals involved in the Magnus series and then doing 
the exponentiation" [49]. Later on, by following a similar approach, Devries 
[62] designed a numerical procedure for determining a fourth order approx- 
imation to the propagator employed in the integration of the single channel 
Schrodinger equation, but it was in the pioneering work [120] where Iserles and 
N0rsett carried out the first systematic study of the Magnus expansion with 
the aim of constructing numerical integration algorithms for linear problems. 
To design the new integrators, the explicit time dependency of each term Q^. 
had to be analyzed, in particular its order of approximation in time to the 
exact solution. 

Generally speaking, the process of rendering the Magnus expansion a prac- 
tical numerical integration algorithm involves three steps. First, the Q series 
is truncated at an appropriate order. Second, the multivariate integrals in 
the truncated series — Vti are replaced by conveniently chosen ap- 
proximations. Third, the exponential of the matrix has to be computed. 
We now briefly consider the first two issues, whereas the general problem of 
evaluating the matrix exponential will be treated in section 5.6. 

We have shown in subsection 2.6 that the Magnus expansion is time symmet- 
ric. As a consequence of cq. (66), if A{t) is analytic and one evaluates its Taylor 
series centered around the midpoint of a particular subintcrval [t„, t„ + /i], then 
each term in Vt^ is an odd function of /i, and thus ^2s+i = 0{h'^'^^^) for s > 1. 
Equivalents, QP^-^l = Q + 0{h^'+^) and fiP^-^l = O + 0{h'^'+^). In other 
words, for achieving an integration method of order 2s (s > 1) only terms up 
to ^2s-2 in the Q series are required [22,120]. For this reason, in general, only 
even order methods are considered. 

Once the series expansion is truncated up to an appropriate order, the mul- 
tidimensional integrals involved have to be computed or at least conveniently 
approximated. Although at first glance this seems to be a quite difficult en- 
terprise, it turns out that their very structure allows one to approximate all 
the multivariate integrals appearing in fl just by evaluating A{t) at the nodes 
of a univariate quadrature [120]. 
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To illustrate how this task can be accomplished, let us expand the matrix A{t) 
around the midpoint, ti/2 = tn + h/2, of the subinterval i„+i], 

°^ / \i 1 d^A(t) 

^(t)=ga,(t-tv2) , where % = JT^U,,, (241) 



1/2 



and insert the series (241) into the recurrence defining the Magnus expansion 
(49)-(50). In this way one gets explicitly the expression of flk up to order 
as 

12 oU 

= (^["0' "0' - «o, ai]) + 0{h') (242) 

^^4 = ^^^[ao,ao,ao,ai] +0(/i'^), 

whereas ^5 = 0{h'^), = 0{h^) and Vt-j = 0{h^). Here we write for clar- 
ity [oii, ai2, . . . , aj(_i, fljj = [flii, [flia, [. . . , [aj,_^, aij . . .]]]. Notice that, as antici- 
pated, only odd powers of h appear in and, in particular, fl2i+i = 0{h'^'^~^^) 
for i > 1. 



Let us denote «j = /i*aj_i. Then [ctj^, ajj, • • • , ctii^i, '^ii] element of order 

/^nH jj^ fact, the matrices {ai}i>i can be considered as the generators (with 
grade i) of a graded free Lie algebra C{ai, a2, ■ ■ ■) [182]. It turns out that it is 
possible to build methods of order p = 2s by considering only terms involving 
cci, . . . , tts in r^. Then, these terms can be approximated by appropriate linear 
combinations of the matrix A{t) evaluated at different points. In particular, 
up to order four we have to approximate 

n^a,-^[a,,a2]+Oih'), (243) 

whereas up to order six the relevant expression is 



11 1 1 

= + Y^«3 - y^[«i,Q;2] + —[a2,a3] + — [ai, ai, ag] (244) 

"210^*^^' ^'^^ 7^^*^^' ^(^^)- 

In order to present methods which can be easily adapted for different quadra- 
ture rules we introduce the averaged (or generalized momentum) matrices 

A»(/i) = i (t - ti/2y A{t)dt = i J^^'^ fA{t + h,2)dt (245) 
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for i — 0, . . . ,s — l.Ii their exact evaluation is not possible or is computation- 
ally expensive, a numerical quadrature may be used instead. Suppose that 6j, 
Ci, {i = 1, . . . ,k), are the weights and nodes of a particular quadrature rule, 
say X (we will use X — G ior Gauss-Legendre quadratures and X — NC for 
Newton-Cotes quadratures) of order p, respectively [2], 



= r^^' A{t)dt = hJ2biAi + 0{hP+'), 



with Ai = A(tn + Cih). Then it is possible to approximate all the integrals A^'^^ 
(up to the required order) just by using only the evaluations A^ at the nodes 
Ci of the quadrature rule required to compute ^4*^°' . Specifically, 



A^'^^hi:b,(cj-iyA,. 



i = 0, . . . , s — 1, 



(246) 



or equivalently, A« = ^E5=i {Qx'^^)-- with {Qx^^).. = bj (9 - |)\ 

In particular, if fourth and sixth order Gauss-Legendre quadrature rules are 
considered, then for s = A; = 2 we have [2] 



61-62--, ci---— , C2-- + — , 



to order four, whereas ior s — k — 3, 



h h ^ h 1 
61=03 = — , 02 = -, Ci = - 

18 9 2 



1 VTE 



10 ' 2' 2+ 10 ' 



to order six, so that 



a/3 Vs 
12 12 



_5_ 

18 



4 
9 



36 " 



18 



\ 24 







_15 

36 

J_ 

24 y 



(247) 



Furthermore, substituting (241) into (245) we find (neglecting higher order 
terms) 



1 J ^- 1 



1 - (-1)^+^' 



Oij, < i < s — 1. 



(248) 
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If this relation is inverted (to order four, s = 2, and six, s = 3) one has 



^(2) = (j,(2)^-l ^ 




9 
4 





-15 





12 





-15 





180 



(249) 



respectively, so that the corresponding expression of ctj in terms of A^*) or Aj 
is given by 



(250) 



Thus, by virtue of (250) we can write Q{h) in terms of the univariate integrals 
(245) or in terms of any desired quadrature rule. In this way, one gets the 
final numerical approximations to Q. Fourth and sixth-order methods can 
be obtained by substituting in (243) and (244), respectively. The algorithm 
then provides an approximation for Y{tn+i) starting from ~ Y{tn), with 

tn+l — tn-\- h. 

The observant reader surely has noticed that, up to order h^, there are more 
terms involved in (242) than those considered in (243) or (244): specifically, 
^CKs in (243) and -^a^ and — ^[tti,Q;4] in (244). The reason is that Jl'^l can 
be approximated by A^'^\ A^^\ A^'^^ up to order and then, these omitted 
terms are automatically reproduced when either A^^\ A'^^\ A^'^^ are evaluated 
analytically or are approximated by any symmetric quadrature rule of order 
six or higher. 

Another important issue involved in any approximation based on the Magnus 
expansion is the number of commutators appearing in Q. As it is already evi- 
dent from (242), this number rapidly increases with the order, and so it might 
constitute a major factor in the overall computational cost of the resulting nu- 
merical methods. It is possible, however, to design an optimization procedure 
aimed to reduce this number to a minimum [23]. For instance, a straightfor- 
ward counting of the number of commutators in (244) suggests that it seems 
necessary to compute seven commutators up to order six in h, whereas the 
general analysis carried out in [23] shows that this can be done with only three 
commutators. More specifically, the scheme 



Ci = [ci;i,o;2], 

C2 = -^[«i,2a3 + Ci] 

Q[6] = + _Lq;3 + _L[_20q;i - as + Ci, as + C2], 



(251) 
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verifies that fi'^l = + 0{h^). Three is in fact the minimum number of com- 
mutators required to get a sixth-order approximation to VL. 

This technique to reduce the number of commutators is indeed vahd for any 
element in a graded free Lie algebra. It has been used, in particular, to obtain 
approximations to the Baker-Campbell-Hausdorff formula up to a given order 
with the minimum number of commutators [18]. 

As an illustration, next we provide the relevant expressions for integration 
schemes of order 4 and 6, which readily follow from the previous analysis. 

Order 4. Choosing the Gauss-Legendre quadrature rule, one has to evaluate 



= A{tn + - ^)h), = A{tn + + ^)h) (252) 

and thus, taking gg'^^ in (247), i?^^) (249) and substituting in (250) we find 

ai = -{Ai+A2), 0:2 = ^(^2 -A). (253) 



Then, by replacing in (243) we obtain 



Qi^]{h) = l{A, + A2)-h'§[A,,A2] 



(254) 



Alternatively, evaluating A at equispaced points, with k — 3 and ci — 0,C2 — 
1/2,C3 = 1; 61 = 63 = 1/6, 62 = 2/3 (i.e., using the Simpson rule to approxi- 
mate//;+'^A(s)cis), 

Ai=A{tn), A2 = A{tn+^), As = A{tn + h) 
we have instead ai = |(^i + 4^2 + A3), a2 = h{A3 — Ai), and then 



QW(/i) = g (Ai + 4A2 + As) -T^[Ai + 4A2 + A3, A3 - A,]. (255) 

It should be noticed that other possibilities not directly obtainable from the 
previous analysis are equally valid. For instance, one could consider 



nW(/i) = -{A, + AA2 + A3) - -[Ai, As]. (256) 

Although apparently more A evaluations are necessary in (255) and (256), this 
is not the case actually, since A^ can be reused at the next integration step. 
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Order 6. In terms of Gauss-Legendre collocation points one has 



Ai^A(tn + {l-^)h), A2^A(tn + lh), As^A(tn + {l + ^)h) 

and similarly we obtain 

ai^hA2, q;2 = — ^(As-Ai), as^ —{A3-2A2 + A1), (257) 

which are then inserted in (251) to get the approximation l^+i = exp(f2[^l)F„. 

If the matrix A{t) is only known at equispaced points, we can use the Newton- 
Cotes (NC) quadrature values with s = 3 and k = 5, bi = = 7/90,62 = 
64 = 32/90,63 = 12/90 and 9 = (j - l)/4, j = 1,...,5. Then, using the 
corresponding matrix Q^^'q from (246) we get 



«i = - '^(^1 + + 28(^2 + A4) + I8A3 
60 

1 

15 

1 

3 

Both schemes involve the minimum number of commutators (three) and re- 
quire three or four evaluations of the matrix A{t) per integration step (observe 
that A5 can be reused in the next step in the Newton-Cotes implementation 
because Ci — and Cg = 1). 

Higher orders can be treated in a similar way. For instance, an 8th-order 
Magnus method can be obtained with only six commutators [23]. Also variable 
step size techniques can be easily implemented [22,118]. 



a2 = ^ (7(^5 - Ai) + 16(^4 - A2)) (258) 

as = \(j{A, + A,) - 4(^2 + A4) - 6A3 



5.4- i From Magnus to Fer and Cayley methods 

For arbitrary matrix Lie groups it is feasible to design numerical methods 
based also on the Fer and Wilcox expansions, whereas for the J-orthogonal 
group (eq. (21)) the Cayley transform also maps the Lie algebra onto the Lie 
group [197] and thus it allows us to build a new class of Lie-group methods. 
Here we briefly show how these integration methods can be easily constructed 
from the previous schemes based on Magnus. In other words, if the solution 
of (223) in a neighborhood of to is written as 
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Y{to + h)^ e^^^^ Yo Magnus (259) 

^QF,{h)^F,{h) ___Yo Fer (260) 

^ QSi(h) g52(ft) . . . g52(ft) gSi{ft) Symmetric Fer (261) 

-(/-^C(/i)) ' + lc(/i)) Fo Cayley (262) 

one may express the functions Fi{h), Si{h) and C{h) in terms of the succes- 
sive approximations to Q and, by using the same techniques as in the previous 
section, obtain the new methods. As the schemes based on the Wilcox fac- 
torization are quite similar as the Fer methods, they will not be considered 
here. 



5.4-2 Fer based methods 

To obtain integration methods based on the Fer factorization (260) one applies 
the Baker-Campbell-Hausdorff (BCH) formula after equating to the Magnus 
expansion (259). More specifically, in the domain of convergence of expansions 
(259) and (260) we can write 

^a{h) ^ gFi(fe) ^F2(h) . . . ^ 

where Fi — Qi is the first term in the Magnus series, F2 — 0{h^) and 
F3 = 0{hJ). Then, a p-th order algorithm with 3 < p < 6, based on the 
Fer expansion requires to compute Fg^^ such that 

Y{tn + h)= e^i(^) e^2''(^) F(t„) + 0{hP+^). (263) 
Taking into account that 

goww ^gF,w e^t'W+0(/,P+i)^ 

we have (Fi = fli) 

Then, by using the BCH formula and simple algebra to remove higher order 
terms we obtain to order four 

1 1 

fJ'^' = -Y^(["i,«2] - - [ai,ai,a2]), (264) 

so that two commutators are needed in this case. A sixth-order method can 
be similarly obtained with four commutators [23]. These methods arc slightly 
more expensive than their Magnus counterpart and they do not preserve the 
time-symmetry of the exact solution. This can be fixed by the self-adjoint ver- 
sion of the Fer factorization in the form (261) proposed in [244] and presented 
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in a more efficient way in [23]. The schemes based on (261) up to order six are 
given by 

Y{tn + h) = e^i('^) e^2'''(/i) e^iW (265) 
with Si — fli/2. A fourth-order method is given by 

^?k^) = -^K,«2] (266) 

(267) 

and a sixth-order one by 



[tti,Q;2j 

[ai, -4a3 + 3si] 

1 



240 



20ai - as + Sl , a2 + ri 



(268) 



To complete the formulation of the scheme, the a, have to be expressed in 
terms of the matrices evaluated at the quadrature points (e.g., equations 
(257) or (258)). 



5.4-3 Cayley-transform methods 

We have seen in subsection 1.2 that for the J-orthogonal group Oj(n) the 
Cayley transform (23) provides a useful alternative to the exponential mapping 
relating the Lie algebra to the Lie group. This fact is particularly important 
for numerical methods where the evaluation of the exponential matrix is the 
most computation-intensive part of the algorithm. 

If the solution of eq. (223) is written as 

Yit) = (/ - ^C(t)) ' (/ + ^C(t)) Yo (269) 
then C{t) e Oj{n) satisfies the equation [113] 

C ^A-^[C,A]- \CAC, t > to, C{to) = 0. (270) 

Time-symmetric methods of order 4 and 6 have been obtained based on the 
Cayley transform (269) by expanding the solution of (270) in a recursive man- 
ner and constructing quadrature formulae for the multivariate integrals that 
appear in the procedure [65,113,163]. It turns out that efficient Cayley based 
methods can be built directly from Magnus based integrators [23]. In partic- 
ular, we get: 
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Order 4: 

CW = - ^in^'^r) = ai - - ^af + Oih% (271) 

where CW = C(/i) + ©(/i^). 
Order 6: 

C[61 = (^/ _ ^(i7[^l)'(/ - ^(^'")')) = C(/i) + 0(/i^). (272) 

Three matrix-matrix products are required in addition to the three commu- 
tators involved in the computation of for a total of nine matrix-matrix 
products per step. 



5.5 Numerical illustration: the Rosen-Zener model revisited 



Next we apply the previous numerical schemes to the integration of the dif- 
ferential equation governing the evolution of a particular quantum two-level 
system. Our purpose here is to illustrate the main issues involved and compare 
the different approximations obtained with both the analytical treatment done 
in section 4 and the exact result. Specifically, we consider the Rosen-Zener 
model in the Interaction Picture already analyzed in subsection 4.1.2. In this 
case the equation to be solved is U[ = Hi(t)Ui, or equivalently, equation (223) 
with Y{t) = Uiit) and coefficient matrix {h= \) 

A{t) = Hi{t) = -iV{s) ((71 cos(^s) - (72 sin(^s)) = -i h{s) ■ a. (273) 

Here V^(s) = Vq/ cosh(s), ^ = cjT and s = t/T. Given the initial condition 
1+) = (1,0)^ at t — — oo, our purpose is to get an approximate value for 
the transition probability to the state |— ) = (0,1)^ ai t — +oo. Its exact 
expression is collected in the first line of eq. (204), which we reproduce here 
for reader's convenience: 

/'„H(f,).(+oo.-oo)P^-Jg^, (274) 

with 7 = nVoT. 

To obtain in practice a numerical approximation to P^x we have to integrate 
the equation in a sufficiently large time interval. We take the initial condition 
at So — ~25 and the numerical integration is carried out until Sf — 25. Then, 
we determine {Ui)i2{sf,So). 

As a first numerical test we take a fixed (and relatively large) time step h such 
that the whole numerical integration in the time interval s e [so, s/] is carried 
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out with 50 evaluations of the vector b(s) for all methods. In this way their 
computational cost is similar. 



To illustrate the qualitative behaviour of Magnus integrators in comparison 
with standard Runge-Kutta schemes, the following methods are considered: 

• Explicit first-order Euler (El): Yn+i = Yn + hA{tn)Yn with tn+i =tn + h 
and h — 1 (solid lines with squares in the figures). 

• Explicit fourth-order Runge-Kutta (RK4), i.e., scheme (227) with h — 
2, since only two evaluations of b(s) per step are required in the linear case 
(soUd Unes with triangles). 

• Second-order Magnus (M2): we consider the midpoint rule (one evalua- 
tion per step) to approximate fli taking h = 1 (dashed lines), i.e., 

y„+i = exp ( - ih b„ • <t) y„ = ^cos(/i6„) I - 5^ . y^. (275) 

with hn = h{tn + h/2) and 6„ = ||b„||. The trapezoidal rule is equally valid 
by considering b„ = (b(t„) + b(t„ + h))/2. 

• Fourth-order Magnus (M4). Using the fourth-order Gauss-Legendre rule 
to approximate the integrals and taking h — 2 one has the scheme (254) 
which for this problem reads 

bi = b(t„ + ci/i), b2 = b(t„ + C2/i), 

d = -(bi + b2) - /i'^i(bi X b2) (276) 
= exp ( - i/i d • <t) y„. 

with ci = |- ^, C2 = |-F^ (dotted lines). 

We choose ^ = 0.3 and ^ = 1, and each numerical integration is carried out 
for different values of 7 in the range 7 e [0,27r]. We plot the corresponding 
approximations to the transition probability in a similar way as in Figure 7 for 
the analytical treatment. Thus we obtain the plots of Figure 14. As expected, 
the performance of the methods deteriorates as 7 increases. Notice also that 
the qualitative behaviour of the different numerical schemes is quite similar 
as that exhibited by the analytical approximations. Euler and Runge-Kutta 
methods do not preserve unitarity and may lead to transition probabilities 
greater than 1 (just like the standard perturbation theory). On the other 
hand, for sufficiently small values of 7 (inside the convergence domain of the 
Magnus series) the fourth-order Magnus method improves the result achieved 
by the second-order, whereas for large values of 7 a higher order method does 
not necessarily lead to a better approximation. 

In Figure 15 we repeat the experiment now taking 7 = 1 and 7 = 2, and for 
different values of ^. In the first case, only the Euler method differs consider- 
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ably from the exact solution and in the second case this happens for both RK 
methods. 




0246 0246 
Y Y 



Fig. 14. Rosen-Zener model: Transition probabilities as a function of 7, with ^ = 0.3 
and ^ = 1. The curves are coded as follows. Solid line represents the exact result; 
El: solid lines with squares; RK4: solid lines with triangles; M2: dashed lines; M4: 
dashed-broken lines (indistinguishable from exact result in left panel). 




Fig. 15. Rosen-Zener model: Transition probabilities as a function of ^, with 7 = 2 
and 7 = 5. Lines are coded as in Figure 14. 

To increase the accuracy, one can always take smaller time steps, but then the 
number of evaluations of A{t) increases, and this may be computationally very 
expensive for some problems due to the complexity of the time dependence 
and/or the size of the matrix. In those cases, it is important to have reliable 
numerical methods providing the desired accuracy as fast as possible or, al- 
ternatively, leading to the best possible accuracy at a given computational 
cost. 
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Fig. 16. Rosen-Zener model: Error in the transition probabilities versus the number 
of evaluations of the Hamiltonian Hj{s) for ^ = 0.3 and 7 = 10 (left panel) and 
7 = 100 (right panel). 

A good perspective of the overall performance of a given numerical integrator is 
provided by the so-called efficiency diagram. This efficiency plot is obtained by 
carrying out the numerical integration with different time steps, corresponding 
to different number of evaluations of A{t). For each run one compares the 
corresponding approximation with the exact solution and plot the error as 
a function of the total number of matrix evaluations. The results are better 
illustrated in a double logarithmic scale. In that case, the slope of the curves 
should correspond, in the limit of very small time steps, to (minus) the order 
of accuracy of the method. 



To illustrate this issue, in Figure 16 we collect the efficiency plots of the 
previous schemes when ^ = 0.3 with 7 = 10 (left) and 7 = 100 (right). We 
have also included the results obtained by several higher order integrators, 
namely the sixth-order RK method (RK6) whose coefficients are collected in 
(229) and the sixth-order Magnus integrator (M6) given by (251) and (257). 
We clearly observe that the Euler method is, by far, the worst choice if accurate 
results are desired. Notice the (double) logarithmic scale of the plots: thus, 
for instance, when 7 = 10 the range goes approximately from 300 to 3000 
evaluations of A{t). Magnus integrators, in addition to providing results in 
SU(2) by construction (as the exact solution), show a better efficiency than 
Runge-Kutta schemes for these examples, and this efficiency increases with 
the value of 7 considered. The implicit Runge-Kutta-Gauss-Legendre methods 
(237) show slightly better results than the explicit RK methods, but still 
considerably worst than Magnus integrators. 

If we compare Magnus integrators of different orders of accuracy, we observe 
that the most efficient scheme is the second order method M2 when relatively 
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(5=0.3,y=10) (^=0.3,Y=1 OO) 




Log (Evaluations) Log (Evaluations) 



Fig. 17. Rosen-Zener model: Same as Figure 16 where we compare the performance 
of the 6th-order Magnus, symmetric-Fer, Cayley and RK6. 

low accuracies are desired. For higher accuracy, however, it is necessary to 
carry out a thorough analysis of the computational cost of the methods for 
a given problem before asserting the convenience of M4 or M6 with respect 
to higher order schemes. For a fixed time step /i, the computational cost of a 
certain family of methods (such as those based on Magnus) usually increases 
with the order. However, if one fixes the number of A{t) evaluations, this 
is not necessarily the case (sometimes higher order methods requires more 
commutators but less exponentials). 

Let us now compare the performance of the Magnus methods with respect to 
other Lie-group solvers, namely Fer and Cayley methods. We repeat the same 
experiments as in Fig. 16 but, for clarity, only the results for the 6th-order 
methods are shown. We consider the symmetric-Fer method given by (265) 
and (268) and the Cayley method (269) with (272) using in both cases the 
Gauss-Legendre quadrature. The results obtained are collected in Figure 17. 
We clearly observe that the relative performance of the Cayley method deteri- 
orates by increasing the value of 7 similarly to the RK6. In spite of preserving 
the qualitative properties, this example shows that for some problems, polyno- 
mial or rational approximations do not perform efficiently. Here, in particular, 
the Magnus scheme is slightly more efficient that the symmetric Fer method, 
although for other problems their performance is quite similar. 



5. 6 Computing the exponential of a matrix 

We have seen that the numerical schemes based on the Magnus expansion 
provide excellent results when applied to equation (223) with coefficient ma- 
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trix (273). In fact, they are even more efficient tlian several Runge-Kutta 
algorithms. Of course, for this particular example the number of A(t) evalua- 
tions is a good indication of the computational cost required by the numerical 
schemes, since the evaluation of exp(O) can be done analytically by means 
of formula (18). In general, however, the matrix exponential has to be also 
approximated numerically and thus the performance of the numerical inte- 
gration algorithms based on the Magnus expansion strongly depends on this 
fact. It may happen that the evaluation of exp(C), where C is a (real or com- 
plex) N X N matrix, represents the major factor in the overall computational 
cost required by this class of algorithms and is probably one of their most 
problematic aspects. 

As a matter of fact, the approximation of the matrix exponential is among 
the oldest and most extensively researched topics in numerical mathematics. 
Although many efficient algorithms have been developed, the problem is still 
far from been solved in general. It seems then reasonable to briefly summarize 
here some of the most widely used procedures in this context. 

Let us begin with two obvious but important remarks, (i) First, one has to 
distinguish whether it is required to evaluate the full matrix exp(C) or only 
the product cxp{C)v for some given vector v. In the later case, special algo- 
rithms can be designed requiring a much reduced computational effort. This 
is specially true when C is large and sparse (as often happens with matrices 
arising from the spatial discretization of partial differential equations), (ii) 
Second, for the numerical integration methods based on ME one has to com- 
pute exp(C(/i)), where C{h) = 0{hP), is a (not too large) step size and 
p > 1. In other words, the matrices to be exponentiated have typically a small 
norm (usually restricted by the convergence bounds of the expansion) . 

In any case, prior to the computation of exp(C), it is significant to have as 
much information about the exponent C as possible. Thus, for instance, if the 
matrix C resides in a Lie algebra, then exp(C) belongs to the corresponding 
Lie group and one has to decide whether this qualitative property has to be 
exactly preserved or constructing a sufficiently accurate approximation (e.g., 
at a higher order than the order of the integrator itself) is enough. Also when 
C can be split into different parts, one may consider a factorization of the 
form exp(C) ~ exp(Ci) exp(C2) • • •exp(C^) if each individual exponential is 
easy to evaluate exactly. 

An important reference in this context is [177] and its updated version [178], 
where up to nineteen (or twenty in [178]) different numerical algorithms for 
computing the exponential of a matrix are carefully reviewed. An extensive 
software package for computing the matrix exponential is Expokit, developed 
by R. Sidje, with Fortran and Matlab versions available [214,215]. In addition 
to computing the matrix-valued function exp(C) for small, dense matrices C, 
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Expokit has functions for computing the vector- valued function exp{C)v for 
both small, dense matrices and large, sparse matrices. 

5.6.1 Scaling and squaring with Fade approximation 

Among one of the least dubious ways of computing exp(C) is by scahng and 
squaring in combination with a diagonal Pade approximation [178]. The pro- 
cedure is based on a fundamental property of the exponential function, namely 

e^ = (e^/^>- 

for any integer j. The idea then is to choose j to be a power of two for 
which exp(C/j) can be rehably and efficiently computed, and then to form 
the matrix (exp(C/j))-' by repeated squaring. If the integer j is chosen as the 
smallest power of two for which ||C||/j < 1, then exp(C / j) can be satisfactorily 
computed by diagonal Pade approximants of order, say, m. This is roughly the 
method used by the built-in function expm in Matlab. 

For the integrators based on the Magnus expansion, as C = 0{h^) with p > 1, 
one usually gets good approximations with relatively small values of j and m. 

As is well known, diagonal Pade approximants map the Lie algebra Oj{n) to 
the J-orthogonal Lie group Oj(n) and thus constitute also a valid alternative 
to the evaluation of the exponential matrix in Magnus-based methods for this 
particular Lie group. More specifically, if S e Oj(n), then t/j2m{tB) e Oj(n) 
for sufficiently small i e R, with 

^2™(A) = (277) 
provided the polynomials Pm are generated according to the recurrence 

Po(A) = l, Pi(A) = 2 + A 
Pm{X) = 2(2m - l)P„,_i(A) + A2p^_2(A). 

Moreover, 7/'2m(A) = exp(A) -|- 0{\'^^^^) and ■02 corresponds to the Cayley 
transform (262), whereas for m = 2, 3 we have 

^.(A) = (l + iA + iA^ + J5A3) / (1 - 1a + iA^ - J5A3) . 

Thus, we can combine the optimized approximations to Q obtained in sub- 
section 5.4 for Magnus based methods with diagonal Pade approximants up 
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to the corresponding order to obtain time-symmetric integration schemes pre- 
serving the algebraic structure of the problem without computing the matrix 
exponential. For instance, the "Magnus-Pade" methods thus obtained involve, 
in addition to one matrix inversion, 3 and 8 matrix-matrix products for order 
4 and 6, respectively. 

Observe that since il^^"! = 0{h) then 

= exp(Q[2"]) + 0{h"'-^'), 

where k = min{m, n}. With m = n we have a method of order 2n. However, 
for some problems this rational approximation to the exponential may be not 
very accurate depending on the eigenvalues of In this case one may take 
m > n, thus giving a better approximation to the exponential and a more 
accurate result by increasing slightly the computational cost of the method. 
Of course, when the norm of the matrix 1]^^"] is not so small, this technique 
can be combined with scaling and squaring 



Instead of using Fade approximants for the exponential of the scaled matrix 
B = C jl^ , Najfeld and Havel [188] propose a rational approximation for the 
matrix function 

R{B) = B coth(B) = B ^5^, (278) 
from which the exponential can be obtained as 

_ H{B) + B 
H{B) - B 

and then iteratively square the result k times to recover the exponential of 
the original matrix C. From the continued fraction expansion of H(B), it is 
possible to compute the first rational approximations as 

7 + ^52 i + ^B^ + ^B^ 

H.iB) = urn = jUb^Ib' 

^ ~ 15 ~ 9 ~ 945 

and so on. Observe that the representation (278) can be regarded as a gener- 
alized Cayley transform of B and thus it also provides approximations in the 
group Oj(n). In [188] the authors report a saving of about 30% in the number 
of matrix multiplications with respect to diagonal Fade approximants when 
an optimal k and a rational approximation for H{B) is used. 



5.6.2 Chebyshev approximation 

Another valid alternative is to use polynomial approximations to the expo- 
nential of C as a whole. Suppose, in particular, that C is a matrix of the 
form C = —irH, with H Hermitian and r > 0, as is the case in Quantum 



107 



Mechanics. In the Chcbyshev approach, the evolution operator exp{—iTH) is 
expanded in a truncated series of Chebyshev polynomials, in analogy with the 
approximation of a scalar function [223]. As is well known, given a function 
F{x) in the interval [—1, 1], the Chebyshev polynomial approximations are op- 
timal, in the sense that the maximum error in the approximation is minimal 
compared to almost all possible polynomial approximations [221]. To apply 
this procedure, one has to previously bound the extreme eigenvalues E^im and 
-Emax of H. Then a truncated Chebyshev expansion of exp(— ix) on the interval 
[rEmin,rEmax] IS Considered: 

m 

exp(-ix) Y,^riPn{x), 

n=0 

where 

^"(^) = i rE -tE ■ j 
with appropriately chosen coefficients c„. Here Tn{x) are the Chebyshev poly- 
nomials on the interval [—1, 1] [2], which can be determined via the recurrence 
relation 

Tn+i{x) = 2xTn{x) - Tn-i{x); Ti{x) = x; To{x) = 1. 
Finally, one uses the approximation 

m 

exp{-zTH) « Yl CnPnirH). (279) 

n=0 

This technique is frequently used in numerical quantum dynamics to compute 
exp{—iTH)ipo over very long times. This can be done with m matrix-vector 
products if the approximation (279) is considered with a sufficiently large 
truncation index m. In fact, the degree m necessary for achieving a specific 
accuracy depends linearly on the step size r and the spectral radius of H [189], 
and thus an increase of the step size reduces the computational work per unit 
step. In a practical implementation, m can be chosen such that the accuracy 
is dominated by the round-off error [145]. This approach has two main draw- 
backs: (i) it is not unitary, and therefore the norm is not conserved (although 
the deviation from unitarity is really small due to its extreme accuracy), and 
(ii) intermediate results are not obtained, since typically r is very large. 



5.6.3 Krylov space methods 

As we have already pointed out, very often what is really required, rather 
than the exponential of the matrix C itself, is the computation of exp(C) 
applied to a vector. In this situation, evaluating e*^ is somehow analogous to 
computing to get the solution of the linear system of equations Cx = b 
for many different 6's: other procedures are clearly far more desirable. The 
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computation of v can be efficiently done with Krylov subspace methods, 
in which approximations to the solution are obtained from the Krylov spaces 
spanned by the vectors {v, Cv, C^v, . . . , C^v} for some j that is typically small 
compared to the dimension of C [71,194]. The Lanczos method for solving 
iteratively symmetric eigenvalue problems is of this form [233]. If, as before, 
we let C = —irH, the symmetric Lanczos process generates recursively an 
orthonormal basis Vm = [f i ■ ■ -Vm] of the mth Krylov subspace Km{H,u) — 
span(M, Hu, . . . , H"^~^u) such that 

where the symmetric tridiagonal mxm matrix Lm = V^HVm is the orthogonal 
projection of H onto Km{H,u). Finally, 

ex.p{-irH)u Vmex.p{-irLm)V^u 

and the matrix exponential exp(— irLj„) can be computed by diagonalizing 

with Dm a diagonal matrix. This iterative process is stopped when 



Pm (exp(-iTL^)) 



< tol 



for a fixed tolerance. Very good approximations are often obtained with rela- 
tively small values of m, and computable error bounds exist for the approx- 
imation. This class of schemes require generally 0{N'^) floating point opera- 
tions in the computation of e'^ v. More details are contained in the references 
[103,105,106]. 



5.6.4 Splitting methods 

Frequently, one has to exponentiate a matrix which can be split into several 
parts which are either solvable or easy to deal with. Let us assume for sim- 
phcity that C — A-\- B, where the computation e'^ is very expensive, but e^ 
and e-^ are cheap and easy to evaluate. In such circumstances, it makes sense 
to approximate e^*^ with e a small parameter, by the following scheme: 

^bl = ^ehmB ^ea^A . . . ^eWB ^ea^A ^ g6(A+B) ^^(gP+l) (280) 

with appropriate parameters flj, 6j. This can be seen as the approximations to 
the solution at t = e of the equation Y' = {A + B)Y by a composition of the 
exact solutions of the equations Y' = AY and Y' = BY at times t = aiS and 
t — biE, respectively. Two instances of this kind of approximations are given 
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by the well known Lie-Trotter formula 




[1] ^ qSA ^eB 



(281) 



and the second order symmetric composition 




[2] ^ g£A/2 g£S g£A/2 



(282) 



referred to as Strang sphtting, Stormer, Verlet and leap-frog, depending on 
the particular area where it is used. 

Splitting methods have been considered in different contexts: in designing sym- 
plectic integrators, for constructing volume-preserving algorithms, in the nu- 
merical integration of partial differential equations, etc. An extensive survey of 
the theory and practice of splitting methods can be found in [27,99,146,165,166] 
and references therein. 

Splitting methods are particularly useful in geometric numerical integration. 

Suppose that the matrix C = A + B resides in a Lie algebra g. Then, obvi- 
ously, exp(C) belongs to the corresponding Lie group Q and one is naturally 
interested in getting approximations also in Q. In this respect, notice that if 
A,B e 0, then the scheme (280) also provides an approximation in Q. It is 
worth noticing that other methods for the approximation of the matrix ex- 
ponential, e.g., Pade approximants and Krylov subspace techniques, are not 
guaranteed to map elements from q to Q. Although diagonal Pade approxi- 
mants map the Lie algebra Oj to the underlying group, it is possible to show 
that the only analytic function that maps sl(n) into the special linear group 
SL(n) approximating the exponential function up to a given order is the expo- 
nential itself [84]. In consequence, diagonal Pade approximants only provide 
results in SL(n) if the computation is accurate to machine precision. 

In [46], Celledoni and Iserles devised a splitting technique for obtaining an 
approximation to exp(C) in the Lie group Q based on a decomposition of 
C e into low-rank matrices Cj e g. Basically, given a n x n matrix C e g, 
they proposed to split it in the form 



k 



such that 



(1) Ci e Q, hi i = 1, . . . ,k. 

(2) Each exp(Ci) is easy to evaluate exactly. 

(3) Products of such exponentials are computationally cheap. 
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For instance, for the Lie algebra 5o{n), the choice 




where ci, . . . , c„ are the columns of C and Cj is the i-th vector of the canonical 
basis of R", satisfies the above requirements (with k — n), whereas in the case 
ol Q — s\{n) other (more involved) alternatives are possible [46]. 

Proceeding as in (281), with 



wc get an order one approximation in e to exp(£C), whereas the symmetric 
composition 



provides an approximation of order two in £, and this can be subsequently 
combined with different techniques for increasing the order. 

With respect to the computational cost, the results reported in [46] show 
that, up to order four in e, this class of splitting schemes are competitive with 
the Matlab built-in function expm when machine accuracy is not required 
in the final approximation. Running expm on randomly generated matrices, 
it is possible to verify that computing exp(C) to machine accuracy requires 
about (20-30)n^ fioating point operations, depending on the eigenvalues of 
C, whereas the 4th-order method constructed from (283) involves (12-15)n^ 
operations when C G 5o(100) [46]. In the case of the approximation of exp(C)'?; 
and V e M'^, the cost of low-rank splitting methods drop down to Kn^ , where 
X is a constant, and thus they are comparable to those achieved by polynomial 
approximants [106]. 

Splitting methods of the above type are by no means the only way to express 
exp(C) as a product of exponentials of elements in q. For instance, the Wei- 
Norman approach (111) can also be implemented in this setting. Suppose that 
dimg = s and let {Xi, X^, . . . , Xg} be a basis of Q. In that case, as we have 
seen (subsection 2.10.3), it is possible to represent exp(tC) for C e g and 
sufficiently small |i| in canonical coordinates of the second kind. 



where the scalar functions arc analytic at t = 0. Although the gi^s are 
implicitly defined, they can be approximated by Taylor series. The cost of 
the procedure can be greatly reduced by choosing adequately the basis and 
exploiting the Lie-algebraic structure [47] . 



— exp(£Ci) exp(£C2) • • • exp(£:Cfc) 



(283) 



QtC _ gPl(t)Xl gS2(t)^2 . . . gi 



111 



Yet another procedure to get approximations of cxp(C) in a Lie-algebraic set- 
ting which has received considerable attention during the last years is based 
on generalized polar decompositions (GPD), an approach introduced in [183] 
and further elaborated in [123,245]. In particular, in [123], by bringing to- 
gether GPD with techniques from numerical linear algebra, new algorithms 
are presented with complexity 0{n^), both when the exponential is applied to 
a vector and to a matrix. This is certainly not competitive with Krylov sub- 
space methods in the first case, but represents at least a 50% improvement on 
the execution time, depending on the Lie algebra considered, in the latter. An- 
other difference with respect to Krylov methods is that the algorithms based 
on generalized polar decompositions approximate the exponential to a given 
order of accuracy and thus they are well suited to exponential approximations 
within numerical integrators for ODEs, since the error is subsumed in that of 
the integration method. For a complete description of the procedure and its 
practical implementation we refer the reader to [123]. 

5.7 Additional numerical examples 

The purpose of subsection 5.5 was to illustrate the main features of the numeri- 
cal schemes based on the Magnus expansion in comparison with other standard 
integrators (such as Runge-Kutta schemes) and other Lie-group methods (e.g.. 
Per and Cayley) on a solvable system. Por larger systems the efficiency analysis 
is more challenging since the (exact or approximate) computation of exponen- 
tial matrices play an important role on the performance of the methods. It 
makes sense then to analyze from this point of view more realistic problems 
where one necessarily has to approximate the exponential in a consistent way. 

As an illustration of this situation we consider next two skew-symmetric ma- 
trices A{t) and Y{0) = I, so that the solution Y{t) of F' = A{t)Y is orthogonal 
for all t. In particular, the upper triangular elements of the matrices A{t) are 
as follows: 



with N = 10. In both cases Y{t) oscillates with time, mainly due to the 
time-dependence of A{t) (in (284)) or the norm of the eigenvalues (in (285)). 

The integration is carried out in the interval t G [0, 10] and the approximate 
solutions are compared with the exact one at the final time t{ = 10 (ob- 
tained with very high accuracy by using a sufficiently small time step). The 
corresponding error at t{ is computed for different values of the time step h. 




l<i< j <N 



(284) 



(285) 
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The Lie-group solvers are implemented with Gauss-Legendre quadratures and 
constant step size. 

First, we plot the accuracy of the different 4-th and 6-th order methods as a 
function of the number of A{t) evaluations. In contrast to the previous exam- 
ples, now there is not a closed formula for the matrix exponentials appearing 
in the Magnus based integrators, so that some alternative procedure must be 
applied. In particular, the computation of e*^ to machine accuracy is done by 
scahng-(diagonal Fade) -squaring, so that the result is correct up to round-off. 
Figure 18 shows the results obtained for the problems (284) and (285) with 
fourth- and sixth-order numerical schemes based on Magnus and Cayley, and 
also explicit Runge-Kutta methods. In the first problem, Magnus and Cayley 
show a very similar performance, which happens to be only slightly better 
than that of RK methods. 

The situation changes drastically, however, for the second problem. Here the 
behaviour of Cayley and RK methods is essentially similar, whereas schemes 
based on Magnus are clearly more efficient. The reason seems to be that Cayley 
and RK methods give poor approximations to the exponential, which, on the 
other hand, has to be accurately approximated, since the eigenvalues of A{t) 
may take large values. 

With respect to symmetric Fer methods, their efficiency is quite similar to 
that of Magnus if the matrix exponentials are evaluated accurately up to 
machine precision. This is so for the matrix (284) even if Fade approximants 
of relatively low order are used to replace the exponentials. 

On the other hand, the efficiency of "Magnus-Fade" methods (we denote by 
MFnm a Magnus method of order n where the exponential is approximated 
by a diagonal Fade of order m, and MPn if n = m) is highly deteriorated 
for the problem (285), although it is always better than the corresponding to 
Cayley schemes. 

To better illustrate all these comments, in Figure 19 we display the error in 
the solution corresponding to (284) and (285) as a function of time in the 
interval t e [0, 100] for h — 1/20 as is obtained by the previous methods. We 
should stress that all schemes require the same number of A evaluations. 

In the right picture the exponentials appearing in the Magnus method are 
computed using a Fade approximant of order six (MF6), of order eight (MF68) 
and to machine accuracy (M6). Observe the great importance of evaluating 
the exponential as accurately as possible for the matrix (285): by increasing 
slightly the computational cost per step in the computation of the matrix 
exponential it is possible to improve dramatically the accuracy of the methods. 
On the contrary, for problem (284) the meaningful fact seems to be that the 
integration scheme provides a solution in the corresponding Lie group. 
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Fig. 18. Efficiency diagram corresponding to the optimized 4-th (circles) and 6-th 
(squares) order Lie-group solvers based on Magnus (solid lines) and Cayley (broken 
lines) , and the standard Runge-Kutta methods (dashed lines) . 




Fig. 19. Error as a function of time (in logarithmic scale) obtained with different 
6-th order integrators for /i = 1/20: (a) problem (284); (b) problem (285). 



5.8 Modified Magnus integrators 



5.8.1 Variations on the basic algorithm 

The examples collected in subsections 5.5 and 5.7 show that the numerical 
methods constructed from the Magnus expansion can be computationally very 
efficient indeed. It is fair to say, however, that this efficiency can be seriously 
affected when dealing with certain types of problems. 
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Suppose, for instance, that one has to numerically integrate a problem defined 
in SU(n). Although Magnus integrators are unconditionally stable in this set- 
ting (since they preserve unitarity up to roundoff^ independently of the step 
size h), in practice only small values of h are used for achieving accurate re- 
sults. Otherwise the convergence of the series is not assured. Of course, the 
use of small time steps may render the algorithm exceedingly costly. 

In other applications the problem depends on several parameters, so that the 
integration has to be carried out for different values of the parameters. In that 
case the overall integration procedure can be computationally very expensive. 

In view of all these difficulties, it is hardly surprising that several modifications 
of the standard algorithm of Magnus integrators had been developed to try to 
minimize these undesirable effects and get especially adapted integrators for 
particular problems. 

One basic tool used time and again in this context is performing a prelimi- 
nary linear transformation, similarly to those introduced in section 2.9. These 
transformations can be carried out either for the whole integration interval 
or at each step in the process. Given an appropriately chosen transformation, 
Fo(^), one factorizes Y{t) = Yiit), where the unknown Yi(t) satisfies the 
equation 

Yl = B{t)Yi (286) 

and B{t) depends on A(t) and Yo(t). This transformation makes sense, of 
course, if ||S(i)|| < ||A(t)|| and thus typically Yq is chosen in such a way 
that the norm of B verifies the above inequality. As a consequence, Magnus 
integrators can be applied on (286) with larger time steps providing also more 
accurate results. 

Alternatively, for problems where in addition to the time step h there is an- 
other parameter (E, say), one may analyze the Magnus expansion in terms 
of h and E. This allows us to identify which terms at each order in the se- 
ries expansion give the main contribution to the error, and design methods 
which include these terms in their formulation. The resulting schemes should 
provide then more accurate results at a moderate computational cost without 
altering the convergence domain. As a general rule, it is always desirable to 
have in advance as much information about the equation and the properties 
of its solution as possible, and then to try to incorporate all this information 
into the algorithm. 

Let us review some useful possibilities in this context. From (241) one has 

A{t) = ao + rai + r^aa H , (287) 

where r = t — ti/2- The first term is exactly solvable (ao = A{ti/2)) and, for 
many problems, it just provides the main contribution to the evolution of the 
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system. In that case it makes sense to take 



YJt) = 6^*"*")"° = e(*~*")'^(*V2) 



and subsequently integrate eq. (286) with 



B{t) = e-(*-*")^(*V2) (A(t) - A(ty2)) e 



(t-t„)A{ti/2) 



This approach has been considered in [115,116], and shows an extraordinary 
improvement when the system is highly oscillatory and the main oscillatory 
part is integrated with Yq. In those cases, the norm of B{t) is considerably 
smaller than ||74(t)||, but B{t) is still highly oscillatory, so that especially 
adapted quadrature rules have to be used in conjunction with the Magnus 
expansion [117,121]. 

In some other problems the contributions from the derivatives can also be 
significant, so that a more appropriate transformation is defined by 



The resulting methods can be considered then as a combination of the Fer 
or Wilcox expansions and the Magnus expansion. This approach has been 
pursued in [61]. 

On the other hand, it is known that several physically relevant systems evolve 
adiabatically or almost-adiabatically. In that case it seems appropriate to con- 
sider the adiabatic picture which instantaneously diagonalizes A{t) (subsection 
2.9). This analysis is carried out in [124,125,151]. In [125] the adiabatic picture 
is used perturbatively, whereas in [124] it is shown that Magnus in the new 
picture leads to significant improvements. 

Alternatively, one can analyze the structure of the leading error terms in order 
to identify the main contribution to the error at each Qi in the Magnus series 
expansion. In most cases they correspond to terms involving only ai and its 
nested commutators with 0:2 Thus, in particular, the standard fourth-order 
method given by (243) can be easily improved by including the dominant 
error term at higher orders, i.e.. 



We recall that using the fourth-order Gauss-Legendre quadrature rule we can 
take ai = ^{Ai+A2), a2 = ^(A2-Ai) with Ai,A2 given in (252). The new 

method requires additional commutators but the accuracy can be improved a 
good deal. This procedure is analyzed in [172], where it is shown how to sum 
up all terms of the form [ai, ai, . . . , ai, 02]. An error analysis in the limit of 
very large values of ||q;i|| is done in [5,160]. 




(288) 



^ai- — [Q!1, Q!2] + Oil, Oil, Oi2] 



1 



[ai, ai, ai, ai, ai, 0:2] + • ■ ■ ■ 



30240 



116 



5.8.2 Commutator- free Magnus integrators 

All numerical methods based on the Magnus expansion appearing in the pre- 
ceding sections require the evaluation of a matrix exponential which contain 
several nested commutators. As we have repeatedly pointed out, computing 
the exponential is frequently the most consuming part of the algorithm. There 
arc problems where the matrix A{t) has a sufficiently simple structure which 
allows one to approximate efficiently the exponential exp{A{ti)), or the ex- 
ponential of a linear combination of the matrix A{t) evaluated at different 
points. In some sense, this is equivalent to have efficient methods to compute 
or to approximate Yq in (288). It may happen, however, that the computation 
of the matrix exponential is a much more involved task due to the presence of 
commutators in the Magnus expansion. For this reason, it makes sense to look 
for approximations to the Magnus expansion which do not involve commu- 
tators whilst still preserving the same qualitative properties. In other words, 
one may be interested in compositions of the form 

*M = exp f p^{s)A{s)ds\ • • • exp f p,{s)A{s)dsj (289) 

where Pi{s) are scalar functions chosen in such a way that \E'["1 = c^(*"+^) -\-0{h^~^^). 
Alternatively, instead of the functions Pi{s), it is possible to find coefficients 
Qij, i = 1, . . . ,m, j = 1, . . . , k such that 

k 

J=l 

is an approximation of the same order. This procedure requires first to com- 
pute Aj = A{tn + Cjh), j = 1, . . . ,k for some quadrature nodes, Cj, of order 
n or higher and, obviously, the coefficients Qij will depend on this choice. 
The process simplifies if one works in the associated graded free Lie algebra 
generated by as in the sequel of eq. (242). Thus, achieving fourth-order 

integrators reduces just to solve the equations for the new coefficients 0^,1, 0^,2 
in 



= exp {a^,! ai + 0^,2 012) • • ■ exp (ai,i ai + ai,i 0:2) (291) 

with the requirement that = exp + 0{h^), where IS given 

by (243). Here the dependence of 0^,1,0^,2 on the coefficients Qij is deter- 
mined through the existing relation between the ctj and the Aj given in 
(250). The order conditions for the coefficients aj,i,at,2 can be easily obtained 
from the Baker-Campbell-Hausdorff formula. As we have already mentioned, 
time-symmetry is an important property to be preserved by the integrators, 
whereas, at the same time, also simplifies the analysis. Scheme (291) is time- 
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symmetric if 



dm+l-i,! — di,!, flm+l-i,2 — — ^1,2, , i — 1,2, ... ,171 (292) 

in which case the order conditions at even order terms are automatically sat- 
isfied. As an illustration, the simple compositions 



= exp Qtti + ^02) cxp Qai - ^^2) (293) 

4^1 = exp(^Q;2) exp(Q;i) exp (^-^as) (294) 

are in fact fourth-order (commutator-free) methods requiring two and three 
exponentials, respectively [25]. In particular, scheme (293), when ai,a2 are 
approximated using the fourth-order Gauss-Legendre quadrature as shown in 
(252) and (253) leads to the scheme 

^J'^ = exp (/i(^2,i^i + ^2,2^2)) exp (h{gi,iAi + ^1,2^2)) (295) 

with ^1^1 = Q2,2 = I + "7I5 Qi,2 = Q2,i = I — Methods closely related 
to the scheme (294) are presented in [10,25,152], where they are apphed to 
the Schrodinger equation with a time-dependent potential. A method quite 
similar to (293) is analyzed in [225] through its application to parabolic ini- 
tial boundary value problems. A detailed study of fourth and sixth order 
commutator- free methods is presented in [28]. 

On the other hand, very often the differential equation (223) can be split into 
two parts, so that one has instead 

Y' ^ [A{t) + B{t))Y, (296) 

where each part can be trivially or very efficiently solved. For instance, the 
Schrodinger equation with a time-dependent potential and, possibly, a time- 
dependent kinetic energy belongs to this class. In principle, the following fam- 
ilies of geometric integrators are specially tailored for this problem: 

I- The commutator- free Magnus integrators (290), which in this case read 

^N = eA^+B^...g^+Bi^ with A = hj2Qi,jA„ B, = hj2Q,,B,. 

(297) 

Assuming that e"^' and e^» are easily computed, then each exponential can 
be approximated by a conveniently chosen splitting method (280) [165] 

gAj+Bj ^ gbs-Bj gds^i . . . gfeiBj gOi^i (298) 
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II- If one takes the time variable in A{t),B(t) as two new coordinates, one 
may use any sphtting method as follows [209]: 

= ^hihB{wi) ^aihA{vi) _ _ _ ^bihB{wi) ^aihA{vi) (299) 

with 

i—l i 

and 60 = 0, A{vi) = A{tn + Vih), B{wi) = B{tn + Wih). 

Both approaches have pros and cons. By applying procedure I we may get 
methods of order 2n with only n evaluations of A(t), B{t) using e.g. Gauss- 
Legendre quadratures, but if m in (297) is large, the number of matrix expo- 
nentials to be computed leads to exceedingly costly methods. The approach 
II, on the other hand, has the advantage of a smaller number of stages, but 
also presents two drawbacks: (i) many evaluations of A{t),B{t) are required 
in general; (ii) for matrices A and B with a particular structure there are 
specially designed splitting methods which arc far more efficient, but these 
schemes are not easily adapted to this situation. 

Next we show how to combine splitting methods with techniques leading to 
commutator free Magnus schemes to design efficient numerical algorithms pos- 
sessing the advantages of approaches I and II, and at the same time general- 
izing the splitting idea (280) to this setting [19,20]. 

The starting point is similar as in previous schemes, i.e. we consider a com- 
position of the form 

=e^' e^'---e^^ (300) 
where the matrices A^ and Bi are 

A = h^PijAj, Bi = h^aijBj, (301) 

with appropriately chosen real parameters pij, aij depending on the coefficients 
of the chosen quadrature rule. Notice that c"^' can be seen as the solution of 
the initial value problem Y' — AiY , Y{tn) = / at tn+i, where A^ — hAi. Of 
course, the same considerations apply to 

In many cases it is convenient to write the coefficients Pij,(Jij explicitly in 
terms of the coefficients q. Following [19] they can be written as 

1=1 ■' 1=1 ■' 

where the coefficients for the matrices R^^\ s — 2,3 are given in (249) and 
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for Qx (whose elements depend on the coefficients bi, Ci for the quadrature 
rule) as shown in (246). 

In this way, the coefficients aij and bij are independent of the quadrature 
choice and can be obtained by solving some order conditions (see [19] for 
more details). 

This procedure allows us to analyse separately particular cases for the matri- 
ces A, B in order to build efficient methods. For instance, in [19] the following 
particular cases are considered: (i) when the matrices A(t), B(t) have a general 
structure; (ii) when they satisfy the additional constraint [B{ti). [B{f j), [B{tk), A{ti)\^ — 
as it happens, for instance, if A corresponds to the kinetic energy and B to 
the potential energy (both in classical or quantum mechanics). 

As an illustration, we consider the following 4th-order 6-stage BAB composi- 
tion 

= c-^« ■ ■ ■ c^^ . (303) 

In Table 1 we collect the coefficients aij,bij to be used in (302) to obtain the 
coefficients Pij,aij to be used in the scheme (303) for two methods, denoted 
by GS6-4 in the general case (whose coefficients aji,&a correspond to Sq in 
[27]) and MN6-4 when [B{ti), [B{tj), [B{tk), A{ti)]]\ = (the coefficients an, b^ 
correspond to SRKN^ in [27]). 

Finally, one has to write the scheme in terms of the matrices Ai,Bi. For in- 
stance, the composition (303) with the 4th-order Gauss-Legendre quadrature 
(i.e. taking Q^^'^^ in (247) and R^'^^ in (249) to obtain the coefficients Pij,(Tij 
in (302)) gives 



A-i = Qoii - V3ai2^ hAi + + V3ai2^ hA2 

Bi = l^^ba - V3bi2j hBi + (^^ba + Vs^^s) hB2. 



(304) 



5. 9 Magnus integrators for nonlinear differential equations 



The success of Magnus methods applied to the numerical integration of linear 
systems has motivated several attempts to adapt the schemes for solving time 
dependent nonlinear differential equations. For completeness we present some 
recently proposed generalizations of Magnus integrators. We consider two dif- 
ferent problems: (i) a nonlinear matrix equation defined in a Lie group, and 
(ii) a general nonlinear equation to which the techniques of section 3.4 can be 
applied. 
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Table 1 

Splitting methods of order 4 for separable non-autonomous systems. GS6-4 is 
intended for general separable problems, whereas MN6-4 can be applied when 
[B{ti),[B{tj),[B{tk),A{ti)]]] = 0. All the coefficients are given in terms of 
ail, 621, a2i, 631 for each method. 



GS 


3-4 


MN6-4 


611 


= 0.0792036964311957 


611 


= 0.0829844064174052 


an 


= 0.209515106613362 


an 


= 0.245298957184271 


621 


= 0.353172906049774 


&21 


= 0.396309801498368 


a2i 


= -0.143851773179818 


021 


= 0.604872665711080 


hi 


= -0.0420650803577195 


631 


= -0.0390563049223486 





= 1/2- (an +a2i) 


hi = 


1 - 2(611 


+ 621 


+ 631) 


an 


= 031 


hi = 


631 






051 


= 021 


hi = 


621 






06I 


= On 


hi = 


611 






012 


= (2aii -1- 2a2i + 031 


-26n 


- 262i)/c 


612 


= (2aii + 2a2i - 26ii 


022 


= 






622 


= (-2an + 6n)/d 


032 


= -an/c 






632 


= 642 = 


a42 


= -032 






652 


= -632 


052 


= -022 






662 


= -622 


062 


= -ai2 






672 


= -612 



c = 12(aii -I- 2a2i + 031 — 26ii -|- 2aii6ii — 2621 + 2aii62i) 
d = 12(2o2i - 611 -I- 2aii6ii - 2a2i6ii - 621 + 2aii62i) 

5.9.1 Nonlinear matrix equations in Lie groups 

As we have already mentioned, the strategy adopted by most Lie-group meth- 
ods for solving the nonlinear matrix differential equation (139), 

r^A{t,Y)Y, Y{0)^Yoeg 

defined in a Lie group Q, whilst preserving its Lie group structure, is to lift 
Y{t) from Q to the underlying Lie algebra q (usually with the exponential 
map), then formulate and numerically solve there an associated differential 
equation and finally map the solution back to ^. In this way the discretization 
procedure works in a linear space rather than in the Lie group. In particular, 
the idea of the so-called Runge-Kutta-Munthe-Kaas class of schemes is to 
approximate the solution of the associated differential equation in the Lie 
algebra g by means of a classical Runge-Kutta method [119,180,181]. 
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To generalize Magnus integrators when A — A{t,Y), an important difference 
witli respect to the Unear case is that now multivariate integrals depend also 
on the value of the (unknown) variable Y at quadrature points. This leads 
to implicit methods and nonlinear algebraic equations in every step of the 
integration [243], which in general cannot compete in efficiency with other 
classes of geometric integrators such as splitting and composition methods. 

An obvious alternative is just to replace the integrals appearing in the non- 
linear Magnus expansion developed in section 3.3 by affordable quadratures, 
depending on the particular problem. If, for instance, we use Euler's method 
to approximate the first term in (142), n^^^{h) = hA{0,Yo) + Oih"^) and O^] jg 
discretized with the midpoint rule, we get the second order scheme 

V2 = hA(^^, et^(0'^'') Yo^ = Q^l (h) + 0{h') 

Yi = e""^ Fo- (305) 

The same procedure can be carried out at higher orders, discretizing consis- 
tently the integrals appearing in Q["*l(/i) for m > 2 [44]. 



5.9.2 The general nonlinear problem 

In principle, it is possible to adapt all methods built for linear problems to 
the general nonlinear non-autonomous equation (159) 

x' = f(i,x), 

or equivalently, the operator differential equation (160), 

As we pointed out in section 3.4, there are two problematic aspects when de- 
signing practical numerical schemes based on Magnus expansion in the non- 
linear case. The first one is how to compute or approximate the truncated 
Magnus expansion (or its action on the initial conditions). The second one is 
how to evaluate the required Lie transforms. For example, to compute the Lie 
transform exp(tLf(y)) acting on y is equivalent to solve the autonomous dif- 
ferential equation x' = f(x) at t — h with x(0) = y, or x(t) = exp(tLf(xo))xo 
where Xq = y can be considered as a set of coordinates. 

Very often, the presence of Lie brackets in the exponent leads to fundamen- 
tal difficulties, since the resulting vector fields usually have very complicate 
structures. Sometimes, however, this problem can be circumvented by using 
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the same techniques leading to commutator-free Magnus integrators in the 
Unear case. In any case, one should bear in mind that the action of the ex- 
ponentials in the methods designed for the linear case has to be replaced by 
their corresponding maps. Alternatively, if the method is formulated in terms 
of Lie transforms, the order of the exponentials has to be reversed, according 
to equation (155). 



Next we illustrate how to numerically solve the problem 

X' = fi(i,x)+f2(i,x) 



(306) 



using the scheme (300) with (304) and the coefficients ajj, bij taken from MN64 
in Table 1. 



(307) 



Let us consider the Duffing equation 

q" + eq' + - 5cos{u}t) 
which can be obtained from the time-dependent Hamiltonian 

H{q,p,t)^T{p,t) + V{q,t)^e-'' ^p^ + e'^ (^^q^ - ^q' - S cos{u;t)q^ (308) 

or equivalcntly from 



d 
— < 
dt 



r{t,p) 

-V'{t,q) 



e 



+ 











e^* {q — q^ + S cos (cut)) 



> . 



(309) 



Notice that this system has already the form (306), each part being exactly 
solvable. In consequence, the splitting method shown in (299) can be used 
here. The procedure is described as Algorithm 1 in Table 2. 



Observe that the leap-frog composition (282) corresponds to m = 2 and 

1 



01 = 02=2' &1 = 1,^2=0. 



(310) 



Since 62 = one stage can be saved (with a trivial modification of the al- 
gorithm) and the scheme is considered as a one stage method. An efficient 
symmetric 5-stage fourth order integrator is given by the coefficients (m = 6) 



bi = 7i- 



(311) 



i = 0, 1, . . . , 6 with 7o = 76 = and 71 = 72 = 74 = 75 = 1/ (4 - 4^3)^ 73 = 
1 - 471. 

Alternatively, we can use the Magnus integrator (303). Since the kinetic energy 
is quadratic in momenta, we can apply the fourth-order method MN6-4. If we 
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Table 2 



Algorithms for the numerical integration of (308) or (309): (Algorithm 1) with 
scheme (299), and (Algorithm 2) with scheme (303). 



Algorithm 1: Standard split 


Algorithm 1: Magnus split 




On = oft», V On = vdm^' 


On = a(t-n)'. Vn = vit-n)'. 


do i = Ij k 




T'iv) = T'iU + Cih v)- V'ia) = V'(U + dh a) 


do ? = 1 77? 

1/ -1. , III' 




^ Ti^ 1 — hn^V^ (f^ i\ 


Ho ? = 1 Tfi 


ta=ta + htti 


Vi{q) = aiiV({q) + ■■■ + UikVl{q) 


Qi = qi-i + hbiT'{tb,Pi) 


Ti{p) = PiiT[{p) H h pikTl.{p) 


tb = tb + hbi 


Pi =Pi-i- hVi{qi-i) 


enddo 


qi = qi-i + hfi{pi) 




enddo 



take the fourth-order Gauss- Legendre quadrature rule for the evaluation of 
the time-dependent function then wc can consider (304), where the coefficients 
ajj, hij are given in Table 1. Here, A(t) plays the role of T'(t, p) and B{t) the role 
of V'{t, q) (they are not interchangeable, otherwise the performance seriously 
deteriorates). The computation of one time step is shown as Algorithm 2 in 
Table 2. 

We take e = 1/20, 5 = 1/4, = 1 and initial conditions g(0) = 1.75, p(0) = 0. 
We integrate up to t = 10 tt and measure the average error in phase space in 
terms of the number of force evaluations for different time steps (in logarithmic 
scale). The results are shown in Figure 20. The scheme MN64 has 6 stages per 
step, but only two time-evaluations. For this reason in the figure we have 
considered as the number of evaluations per step both two and six (left and 
right curves connected by an arrow). It is evident the superiority of the new 
splitting Magnus integrators for this problem. If the time-dependent functions 
dominate the cost of the algorithm the superiority is even higher. Surprisingly, 
the method shows better stability than the leap-frog method, which attains 
the highest stability possible among the splitting methods for autonomous 
problems. 



6 Some applications of the numerical integrators based on ME 

In this section we collect several examples where the numerical integration 
methods based on the Magnus expansion have been applied in the recent 



124 



1.5 2 2.5 3 3.5 

LOG(N. EVAL.) 

Fig. 20. Average error versus number of force evaluations in the numerical inte- 
gration of (309) using second and fourth order symplectic integrators for general 
separable systems (S2 corresponds to the second order leapfrog method with coef- 
ficients (310) and SU54 to the fourth order method with coefficients (311)) and the 
fourth order symplectic Runge-Kutta-Nystrom method MN64 with initial condi- 
tions q{0) = 1.75, p{0) = and e = 1/20, 6 = 1/4, u = 1. 

literature. Special attention is dedicated to the numerical integration of the 
Schrodinger equation, since the Magnus series expansion has been extensively 
used in this setting almost since its very formulation. The time-independent 
Schrodinger equation can be considered as a particular example of a Sturm- 
Liouville problem, so we also review the applicability of Magnus based tech- 
niques in this context. Then we consider a particular nonlinear system (the 
differential Riccati equation) which can be, in some sense, linearized, so that at 
the end one may work with finite-dimensional matrices. Finally, we summarize 
a recent but noteworthy application: the design of new classes of numerical 
schemes for the integration of stochastic differential equations. 



6.1 Case study: numerical treatment of the Schrodinger equation 

Before embarking ourselves in the use of numerical methods based on the Mag- 
nus expansion in the integration of the Schrodinger equation, let us establish 
first the theoretical framework which allows one to use numerical integrators 
in this setting for obtaining approximate solutions in time and space. 
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6.1.1 Time- dependent Schrddinger equation 



To keep the treatment as simple as possible, we commence by considering the 
one-dimensional time-dependent Schrodinger equation {h — 1) 

i^^Pit, x) = Htlj{t, x) = -^^^(^, + V{x)i;{t, x), (312) 

with initial condition iIj{0,x) = ipo{x). If we look for a solution of the form 
ip{t,x) = (j){t) <f{x), it is clear that, by substituting into (312), one gets (f){t) = 
Q-ttE^ where E is a. constant and (p{x) is the solution of the second order 
differential equation 

- ^ + V{x)ip = Eip. (313) 

li E > V the solution is oscillatory, whereas if E < V the solution is a 
linear combination of exponentially increasing and decreasing functions. For 
bounded problems this last condition always takes place at the boundaries. 
Since 

J \il;{x,t)\'^dx ^ J \(fi{x)\'^dx < oo, (314) 

it is clear that the exponentially increasing solutions have to be cancelled, and 
this can only occur for certain values of the constant E, which are precisely 
the eigenvalues of the problem. 

Let us assume that the system has only a discrete spectrum and denote by 
{E„, </?n}^0' with Ei < Ej, i < j, the complete set of eigenvalues and associ- 
ated eigenvectors. It is well known that we can take {(pn}'^=Q as an orthonormal 
basis and, since (312) is linear, any solution can be written as 

oo 

V^(t,a:) = ^c„e-*^"(^„(x). (315) 

n=0 

Using the standard notation for the inner product, one has 

{iPnix)m, x)) = I iplix) i^it, x)dx = c„ e-^*^" (316) 

and 

is the probability to find the system in the eigenstate so that knP — 1- 
The energy is given by 

„ oo oo 

8 = {iP\H\,p) = / r{t,x)H,p{t,x)dx =Y.Y. <C^^n,m, (317) 

n=l m=l 

where 

Hn,m = {m\H\n) = {(pm\H\ipn). 
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In general, the coefficients c„ decrease very fast with n and. in some cases, 
the system allows only a finite number of states. In that situation, one may 
consider the Schrodinger equation as a finite dimensional linear system where 
the Hamiltonian is a matrix with elements i^n.m- This is precisely the case for 
the examples examined in section 4. 

When the Hamiltonian is explicitly time-dependent, this procedure is not 
longer valid. Instead one may use some alternative techniques which we now 
briefly review. 

(i) Spectral Decomposition. Let us assume that the system is perturbed with 
a time-dependent potential, i.e., equation (312) takes the form 

x) = Hmit, x) = {f + V{t)) i{j{t, x), (318) 

where 

In this case we cannot use separation of variables. However, since {(fin} is a 
complete bases we can still write the solution as 

d-l 

.Pit, x) ~ 5] c„(t) e-^*^" <fin{x), (319) 

Tt=0 

where and are the exact eigenvalues and eigenfunctions when V = 0, 
and the complex coefficients c„ give the probability amplitude to find the 
system in the state cpn (Z^n |cn(^)P — 1 for all t). Then, substituting (319) into 
(318) we obtain the matrix equation 

ij^c{t)^H{t)cit), c(0)=co, (320) 

where c = (cq, . . . , Q-i)^ G and H e ([^dxd jg Hermitian matrix 
associated to the Hamiltonian 

{H{t))ij = {(pi\H{t) - Ho\(pj) e^(^'-^^)*, i,j^l,...,d 

and ^0 — H{t = 0). Given the initial wave function -0(0,2;), the components 
of Co are determined by Co,^ = {(fii\ip{0,x)). 

Obviously, any complete basis can be used in this case, although the norm 
of the matrix H{t) may depend on the choice. In addition, the number of 
basis elements, i.e. the minimum dimension d necessary to obtain a sufficiently 
accurate result, also depends on the chosen basis. 

(ii) Space discretization. This procedure intends to take profit of the structure 
of the Hamiltonian H in (318): V is diagonal in the coordinate space and T is 
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diagonal in the momentum space. Let us assume that the system is defined in 
the interval x G [a;o,a;y] with periodic boundary conditions. We can then split 
this interval in d parts of length Ax = {xf —XQ)/d and consider c„ = ip{t, Xn) 
where Xn = xq + nAx, n = 0,1, . . . ,d — 1. Then a finite dimensional linear 
equation similar (but with a different coefficient matrix H) to equation (320) 
results. Since V is diagonal in the coordinate space and T is diagonal in 
momentum space, it is possible to use complex Fast Fourier Transforms (FFTs) 
for evaluating the products He, where Tipit, Xn) = J-~^DTJ-'ip{t, Xn), and Dt 
is a diagonal operator. 

We thus see that whatever the procedure used (spectral decomposition or 
space discretization) , one ends up with a linear equation of the form 

t^{t)^H{t)m, m-^o (321) 

where now iplt) represents a complex vector with d components which approx- 
imates the (continuous) wave function. The computational Hamiltonian H{t) 
appearing in (321) is thus a space discretization (or other finite-dimensional 
model) of H{t) — T + V{t). Numerical difficulties come mainly from the un- 
bounded nature of the Hamiltonian and the highly oscillatory behaviour of 
the wave function. 

It is at this point when numerical algorithms based on the Magnus expansion, 
such as they have been formulated in previous sections, come into play for 
integrating in time the linear system (321). To put them in perspective, let us 
introduce first some other numerical methods also used in this context. Our 
exposition is largely based on the reference [156]. 



6.1.1.1 The implicit midpoint rule The approximation to the solution 
of (321) provided by this scheme is implicitly defined by 

^tn+^^ = i/(tn+l/2) ^(^n+1 + A), (322) 

where tn+1/2 = |(^n+i+^n)- Here and in the sequel, for clarity, we have denoted 
by At the time step size and tn — nAt. Alternatively, 

ipn+i^r{-iAtH {tn+1/2)) ipn, with r{z)^\^^. (323) 

Observe that, as r is nothing but the Cayley transform, the numerical prop- 
agator is unitary and consequently the Euclidean norm of the discrete wave 
function is preserved along the evolution: = \\ipn\\- This is a crucial 

qualitative feature the method shares with the exact solution, contrarily to 
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other standard numerical integrators, such as exphcit Runge-Kutta methods. 
From a purely numerical point of view, the algorithm is stable for any step 
size At. 

Another useful property of this numerical scheme is time-symmetry: exchang- 
ing in (323) n by n + 1 and At by —At we get the same numerical method 
again. Equivalent ly, r{—z) = r{z)~^, exactly as the exponential e^. 

With respect to accuracy, it is not difficult to show that, if H{t) is bounded 
and sufficiently smooth, the error verifies 



uniformly for nAt in a time interval [O.t/]. In other words, the implicit mid- 
point rule is a second-order method. It happens, however, that the constant 
in the term 0{At^) depends on bounds of H' and H" and on the norm of 
the third derivative of the solution ip. Since, in general, the wave function is 
highly oscillatory in time, this time derivative can become large, and so the 
use of very small time steps is mandatory. 

6.1.1.2 The exponential midpoint rule Another possibility to get ap- 
proximate solutions of (321) consists in replacing r{z) by exp(2;) in (323): 



Now, instead of solving systems of linear equations as previously, one has to 
compute the exponential of a large matrix times a vector at each integration 
step. In this respect, the techniques reviewed in subsection 5.6 can be efficiently 
implemented. The exponential midpoint rule (325) also provides a unitary 
propagator and it is time-symmetric. In addition, the error satisfies the same 
condition (324), but now the constant in the 0{At'^) term is independent 
of the time derivatives of ip under certain assumptions on the commutator 
[H{t),H{s)] [104]. As a consequence, much larger time steps can be taken to 
achieve the same accuracy as with the implicit midpoint rule. 

6.1.1.3 Integrators based on the Magnus expansion The method 
(325) is a particular instance of a second order Magnus method when the 
integral J^^ H{s)ds is replaced by the midpoint quadrature rule. In fact, we 
have already used it in (275). Obviously, if higher order approximations are 
considered, the accuracy can be enhanced a great deal. This claim has to 
be conveniently justified, however, since the order of the numerical methods 
based on Magnus has been deduced when ||AtiJ(t)|| and is obtained 
by studying the remainder of the truncated Magnus series. In the Schrodinger 
equation, on the other hand, one has to cope with discretizations of unbounded 



Un-1p{Q\\=0{At^) 



(324) 



Ipn+l = exp{-iAtH{tn+i/2)) Ipn- 



(325) 
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operators, so in principle it is not evident how the previous resuhs on the order 
of accuracy apply in this context. In [104], Hochbruck and Lubich analyse in 
detail the application of the fourth-order Magnus integrator (254) to equation 
(321), showing that it works extremely well even with step sizes for which the 
corresponding \\AtH(t) \\ is large. In particular, the scheme retains fourth order 
of accuracy in At independently of the norm of H{t) when H{t) = T + V(t), T 
is a discretization of (with maximum eigenvalue -Emax ~ (Aa;)~^) and 

V{t) is sufficiently smooth under the time step restriction /S.t^/E^^ < Const. 
This is so even when there is no guarantee that the Magnus series converges 
at all. 



6.1.1.4 Symplectic perspective The evolution operator corresponding 
to (318) is not only unitary, but also symplectic with canonical coordinates 
and momenta Re('^) and Im('^), respectively. If we carry out a discretization 
in space, this symplectic structure is inherited by the corresponding equation 
(320). It makes sense, then, to write c = q -|- ip and consider the equations 
satisfied by q, p e M'^, namely 



H(t)p, p' = -H(i)q, 



(326) 



which can be interpreted as the canonical equations corresponding to the 
Hamiltonian [95] 

n{t, q, p) = p^H(t)p + q^H(t)q. (327) 

Denoting z = (q, p)^, it is clear that 

z' = {Ait) + B(t)) z 



where 



A{t) 




(328) 



For this system it is possible, therefore, to apply the commutator-free Magnus 
integrators constructed in subsection 5.8.2. In addition, one has 

[B,[B,[B,A]]] = [A, [A, [A,B]]]=0, 

and this property allows us to use especially designed and highly efficient 
integration methods [20]. 



6.1.2 Time-independent Schrddinger equation 

Restricting ourselves to the time-independent Schrodinger equation (313), we 
next illustrate how Magnus integrators can in fact be used to compute the 
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discrete eigenvalues defined by the problem. Although only the Schrodinger 
equation in a finite domain is considered, 



- ^ + V{x)lp = Xip, xe {a, h) (329) 

the procedure can be easily adapted to other types of eigenvalue problems, in 
which one has to find both A = and (p. Here it is assumed that the potential 
is smooth, V e C'^{a, h) and, for simphcity, <^(a) = (^(6) = 0. 

Under these assumptions, it is well known that the eigenvalues are real, distinct 
and bounded from below. The problem (329) can be formulated in the special 
linear group SL(2), 



dx 








y, xe(a, 6), where y — {ip,d(p/dx)'^, 



(330) 

so that the Magnus expansion can be applied in a natural way. As usual, rather 
than approximating the fundamental solution of (330) in the entire interval 
(a, b) by exp(Jl), the idea is to partition the interval into A^" small subintervals, 
and then apply a conveniently discretized version of the Magnus expansion. 
In this way, the convergence problem no longer restricts the size {b — a) [172]. 

For the sake of simphcity, let us consider the fourth-order method (254). Writ- 
ing 

K,l = V{Xn + (^ - ^)h), Vn,2 = V{Xn + (^ + ^)h), 

Z b ZD 
where h = {b — a) /N and Xn = a + hn, we form 



(7-n(A) = 



( -^/i2(K„,i-K,2) h 



for n = 0, 1, . . . , A" — 1. Then, the fourth-order approximation to the solution 
of (330) Bit X — b is 

y{b) = e'^^-iW • • • e^i W e^o(^) y(a) (331) 

and the values of A are obtained from (331) by using repeatedly the expression 
of the exponential of a traceless matrix, eq. (24), and requiring that ip{a) = 
ip{b) — 0. The resulting nonlinear equation in A can be solved, for instance, by 
Newton-Raphson iteration, which provides quadratic convergence for starting 
values sufficiently near the solution [172]. 

Although by construction this procedure leads to a global order of approxi- 
mation 0{h^) if a pth-order Magnus method is apphed, it turns out that the 
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error also depends on the magnitude of the eigenvalue. Specifically, the error 
in a pth-order method grows as [172], and thus one expects 

poor approximations for large eigenvalues. This difficulty can be overcome up 
to a point by analyzing the dependence on A of each term in the Magnus 
series and considering partial sums of the terms carrying the most significant 
dependence on A. For instance, it is possible to design a sixth-order Magnus 
integrator for this problem with error (9(/?7A). which therefore behaves like a 
fourth-order method when h^X a; 1, whereas the standard sixth-order Magnus 
scheme, carrying an error of 0{hJX'^), reduces to an order-two method [172]. In 
any case, getting accurate approximations when [A| — > oo is more problematic 
[128]. 



6.2 Sturm-Liouville problems 



The system defined by (329) with boundary conditions ^{a) = ip{h) = is just 
one particular example of a second order Sturm-Liouville problem [199,246]. 
It is thus quite natural to try to apply Magnus integrators to more general 
problems within this class. 

A second order Sturm-Liouville eigenvalue problem has the form 

— (^{x)-^{x)^ + q{x)y{x) = Xr{x)y{x) on (a, 6) (332) 

with separated boundary conditions which commonly have the form 

A,y(a) + A2p(a)y\a) = B,y(b) + B2p(b)y'(b) = (333) 

for given constants Ai, Bi and functions p(a;), q{x) and r{x). Solving this prob- 
lem means, of course, determining the values A„ of A for which eq. (332) has 
a nontrivial (continuously differentiable square integrable) solution y„(x) sat- 
isfying equations (333) [246,7]. 

These and other higher order Sturm-Liouville problems can be recasted as a 
linear matrix system of the form 

Y' = {XB + C{x)) Y (334) 

by transforming to the so-called compound matrix or modified Riccati vari- 
ables [96,128]. Here S is a constant matrix. When generahzing the above treat- 
ment based on the Magnus expansion to this problem, there is one elementary 
but important remark worth to be stated explicitly: unless the differential 
equation (334) has the same large X-asymptotics as some differential equation 
with x-independent coefficients, then it will be impossible to develop a Magnus 
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method which accurately approximate its solutions for large A [128]. The rea- 
son is that a Magnus method approximates the solution by a discrete solution 
calculated using a formula of the form Y{xn+i) = exp((T„(A))F(a;„); in par- 
ticular, on the first step {xo,Xi), the differential equation is approximated by 
one in which the coefficient matrix is replaced by the x-independent matrix 

(7o(A)/(xi - Xo). 



In consequence, the attention should be restricted to systems for which it is 
known that a suitable constant-coefficient system provides the correct asymp- 
totics. This is the case, in particular, for equation (329), and more generally 
for linear equations of order 2n in which the (2n — l)st derivative is zero, such 
as 



2n-2 

i=o 



Here the asymptotics are determined by the equation 1)"?/(2") = Xy [187]. 
Even then, the methods developed in [172] for equation (329) and implemented 
for systems with matrices of general size in [128] require a A-dependent step 
size restriction of the form h < 0(|A|~^/^) in order to be defined. Neverthe- 
less, the analysis carried out in [128] shows that the fourth order Magnus 
integrator based on a two-point Gaussian quadrature appears to offer signifi- 
cant advantages over conventional methods based on power series and library 
routines. 



Magnus integrators have also been successfully applied in the somewhat re- 
lated problem of computing the Evans function for spectral problems arising 
in the analysis of the linear stability of travelling wave solutions to reaction- 
diffusion PDEs [5]. In this setting, Magnus integrators possess some appealing 
features in comparison, for instance, with Rungc-Kutta schemes: (1) they are 
unconditionally stable; (2) their performance is superior in highly oscillatory 
regimes and (3) their step size can be controlled in advance. Items (2) and (3) 
are due to the fact that error bounds for Magnus methods depend only on low 
order derivatives of the coefficient matrix, not (as for Runge-Kutta schemes) 
on derivatives of the solution. Therefore, performance and, correspondingly, 
the choice of optimal step size remain uniform over any bounded region of 
parameter space [5]. 
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6.3 The differential Riccati equation 



Let us consider now the two-point boundary value problem in the t variable 
defined by the linear differential equation 



y = 



Vy2. 



A{t) B{t) 
C{t) D{t) 




0<t<T 



(335) 



with separated boundary conditions 



{Kn Ku) 




7i, {K21 K22) 



t=o 




72- 



(336) 



t=T 



Here A e , B e O^^, C e 0""^ , D E C^*' , whereas yi, 72 G C^, ya, 71 G 
C and the matrices Kij have appropriate dimensions. We next introduce the 
time-dependent change of variables (or picture) y = Fo(^) w, with 



/ 



X{t) I, 







(337) 



and choose the matrix X e C^^? go as to ensure that in the new variables 
w = YQ^{t)y the system assume the partly decoupled structure [63] 



w 





1 








V 



A + BX B 
O D-XB 



1 






I W2 / 



(338) 



together with the corresponding boundary conditions for w. It turns out that 
this is possible if and only if X{t) satisfies the so-called differential Riccati 
equation [64] 

X' = C{t) + D{t)X-XA{t)-XB{t)X, X(0)=Xo (339) 
for some Xq. By requiring 

Xo^-K:[^Kr^, (340) 

then the boundary conditions (336) also decouple as 

(O Xi2)w(0) = 71, (K2i + K22X{T) ir22)w(r) = 72. (341) 

Here we assume without loss of generality that K12 is invertible. In this way, 
the original boundary value problem can be solved as follows [63,64]: (i) solve 
equation (339) with initial condition (340) from t — Q to t — T; (n) solve 
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the W2-equation in (338) and (341), also from zero to T; (iii) solve the wi- 
equation in (338) from t = T to t = and recover y = YQ{t)w. In other 
words, the solution of the original two-point boundary value problem can be 
obtained by solving a sequence of three different initial value problems, one 
of which involves the nonlinear equation (339). Obviously, steps (i) and (iii) 
can be solved using numerical integrators based on the Magnus expansion. It 
could be perhaps more surprising that these algorithms can indeed be used to 
integrate the Riccati equation in step (ii). 

Although the boundary value problem (335) is a convenient way to introduce 
the differential Riccati equation (339), this equation arises in many fields of 
science and engineering, such as linear quadratic optimal control, stability 
theory, stochastic control, differential games, etc. Accordingly, it has received 
considerable attention in the literature, both focused to its theoretical aspects 
[17,200] and its numerical treatment [63,133,211]. 

In order to apply Magnus methods to solve numerically the Riccati equation, 
we first apply the transformation 

X{t) = V{t)W-\t), (342) 

with V e Cf W e C«^« and V(0) = Xq, W{0) = Ig, in the region where 
W{t) is invertible. Then eq. (339) is equivalent to the hnear system 

Y' = S{t) Y{t) , r (0) = I I (343) 



^0 



with 



. (344) 
Vit) ) \c(t) D{t) ) 

SO that the previous Magnus integrators for linear problems can be applied 
here. Apparently, this system is similar to (335), but now we are dealing with 
an initial value problem and y is a matrix instead of a vector. 

When dealing in general with the differential Riccati equation (339), it is 
meaningful to distinguish the following three cases: 

(i) The so-called symmetric Riccati equation, which corresponds to q — p, 
D{t) = —A{t)'^ real, and B{t), C{t) real and symmetric matrices. In this 
case, the solution satisfies X'^ = X. It is straightforward to show that this 
problem is equivalent to the treatment of the generalized time-dependent 
harmonic oscillator, described by the Hamiltonian function 

H = \p'B{t)p + p^A(t)q - ^q^C(t)q. 
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The approximate solution attained by Magnus integrators when applied to 
(343)- (344) can be seen as the exact sohition corresponding to a perturbed 
symplectic matrix S{t) ~ S{t). In other words, we are solving exactly a 
perturbed Hamiltonian system so that the approximate solution, X, will 
shares several properties of the exact solution, in particular — X. 

(ii) The linear non-homogeneous problem 

X' = D(t)X + C{t) (345) 

corresponds to the particular case A — Q and B — Q va. (339). 

(iii) The problem 

X' = D{t)X + XA{t) (346) 

is recovered from (339) by taking C = and i? = 0. It has been treated in 
[112] by developing an ad hoc Magnus- type expansion. Notice that the case 
p = q, D — —A corresponds to the hnear isospectral system (238). 

6.4 Stochastic differential equations 

In recent years the use of stochastic differential equations (SDEs) has become 
widespread in the simulation of random phenomena appearing in physics, 
engineering, economics, etc, such as turbulent diffusion, polymer dynamics 
and investment finance [35] . Although models based on SDEs can offer a more 
realistic representation of the system than ordinary differential equations, the 
design of effective numerical schemes for solving SDEs is, in comparison with 
ODEs, a less developed field of research. This fact notwithstanding, it is true 
that recently new classes of integration methods have been constructed which 
automatically incorporate conservation properties the SDE possesses. Since 
some of the methods are based precisely on the Magnus expansion, we briefly 
review here their main features, and refer the reader to the more advanced 
literature on the subject [35,141,196]. 

A SDE in its general form is usually written as 

d 

dy{t) = go{t,y{t)) dt + Y.9j{t.y{t)) dW,{t), y(0) = |/o, ye M'", 

(347) 

where Qj, (j > 0), are m- vector- valued functions. The function qq is the deter- 
ministic continuous component (called the drift coefficient), the Qj, (j > 1), 
represent the stochastic continuous components (the diffusion coefficients) 
and Wj are d independent Wiener processes. A Wiener process W (also called 
Brownian motion) is a stochastic process [35] satisfying 

W(Q) = 0, E[W{t)] = 0, Yar[W{t) - W{s)] ^t-s, t>s 
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which has independent increments on non- overlapping intervals. In other words, 
a Wiener process is normally distributed with mean or expectation value E 
equal to zero and variance t. 

Equation (347) can be written in integral form as 

y{t)^yo+ fgo{s,yis))ds + Yl f gj{s,y{s))dWj{s). (348) 
JO j=i 

The d integrals in (348) cannot be considered as Riemann-Stieltjes integrals, 
since the sample paths of a Wiener process are not of bounded variation. In 
fact, if different choices are made for the point (in the subintervals 
of a given partition) where the function is evaluated, then the approximating 
sums for each gj, 

N 

T.9jiri,y{Ti)){Wj{U) - WjiU.i)), n = eu + (1 - ^)ti-i, (349) 
1=1 

converge (in the mean-square sense) to different values of the integral, depend- 
ing on the value of [34] . Thus, for instance, 

[ W{t)dW{t) = \{W\h) - W\a)) + {e-\){h- a). 

If ^ = 0, then Ti — ti-i (the left-hand point of each subinterval) and the 
resulting integral is called an Ito integral; if ^ = 1/2 (so that the midpoint is 
used instead) , one has a Stratonovich integral. These are the two main choices 
and, although they are related, the particular election depends ultimately 
on the nature of the process to be modeled [34]. It can be shown that the 
Stratonovich calculus satisfies the Riemann-Stieltjes rules of calculus, and 
thus it is the natural choice here. 

When dealing with numerical methods for solving (347), there are two ways of 
measuring accuracy [35]. The first is strong convergence, essential when the aim 
is to get numerical approximations to the trajectories which are close to the 
exact solution. The second is weak convergence, when only certain moments of 
the solution arc of interest. Thus, if y„ denotes the numerical approximation to 
y{tn) after n steps with constant step size h = (tn — to)/n, then the numerical 
solution y converges strongly to the exact solution y with strong global order 
p if there exist C > (independent of h) and S > such that 

E[\\yn-y{tn)\\]<ChP: he (0,5). 

It is worth noticing that p can be fractional, since the root mean-square order 
of the Wiener process is /i^/^. One of the simplest procedures for solving (347) 
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numerically is the so-called Euler-Maruyama method [164] , 

d 

Vn+l ^Vn + Yl Jj9j(tn, Vn) , (350) 

where 

h = tn+l - tn, Jo = h, Jj = Wj{tn+l) - Wj{tn), j = l,...,d. 

This scheme turns out to be of strong order 1/2. Here the Jj can be computed 
as \/hNj, where the Nj are A^(0, 1) normally distributed independent random 
variables [34]. 

For the general non-autonomous linear Stratonovich problem defined by 

dy = Go{t)y dt + Yl Gjit) V dWj, y{0) ^yoeW^ (351) 

the Magnus expansion for the deterministic case can be extended in a quite 
straightforward way. It is well worth noticing that in equation (351), even 
when the functions Gj are constant, there is no explicit solution [141], unless 
all the Gj, j > 0, commute with one another, in which case it holds that 

y{t) = exp (cot + Y^ GjWjit)j y^. (352) 

In many modeling situations, however, there is no reason to expect that the 
functions Gj associated with the Wiener processes commute. If for simplicity 
we only consider the autonomous case and write 

d 

G{t) = Godt + Y^ Gj dWj{t) , 

then (351) can be expressed as 

dy^G{t)ydt, y(0) = yo 

and thus one can formally apply the Magnus expansion to this equation to 
get y{t) — exp{Q{t))yo. The first term in the series reads in this case 

/* G(s) ds = [*Gods + y t G. dWAs) ^G^t + Y G^J., 
Jo Jo ^^Jo ^ 



where now 



Jj^ fdW(s)^Wj(t)-Wj(0). 

J 
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By inserting these expressions into the recurrence associated with the Magnus 
series, Burrage and Burrage [34] show that 



1 ^ 



nit) = J2 GjJj + ^ E E G,]{J,, - J,,) (353) 

j=0 i=0 j=i+l 

d d d / -j^ s 

+ E E E [^«' t^i' Gk]] ( l^iJkji - Jjki) + TT^MJjk - Jkj) ) H ) 

where the multiple Stratonovich integrals are defined by 

jj^j^-jxt) = f r ■■■ r ^^ji(^i) ■ ■ ■ d^n{^i)^ J^ e {0, 1, . . . , d}. 

t/ 

(354) 

Since not all the Stratonovich integrals are independent, one has to compute 
only d{d + l)(d + 5)/6 stochastic integral evaluations to achieve strong order 
1.5 with the expression (353) [34]. If, on the other hand, Q{t) is truncated 
after the first set of terms, then the resulting numerical approximation 

y{t) = exp ^GjJ^ yo 

has strong order 1/2, but leads to smaller error coefficients than the Euler- 
Maruyama method (350) [34]. Furthermore, the error becomes smaller as the 
GiGj terms get closer to commuting and the scheme preserves the underlying 
structure of the problem. 

One should notice at this point that equation (347) (or, in the linear case, 
equation (351)), has formally the same structure as the nonhnear ODE (171) 
appearing in control theory. Therefore, the formalism developed there to get 
the Chen-Fliess series can be applied here with the alphabet I — {0,1, ... ,d} 
and the integrals 

J^l?j (t) = j\{s)ds, (t) = j\{s)dW,{s), I > 1, 

since the Stratonovich integrals satisfy the integration by parts rule. In other 
words, one can obtain the corresponding Magnus expansion for arbitrary (lin- 
ear or nonlinear) stochastic differential equations simply by following the same 
procedure as for deterministic ODEs. 

With respect to nonlinear Stratonovich stochastic differential equations, it 
should be remarked that the use of Lie algebraic techniques as well as the 
design of Lie group methods for obtaining strong approximations when the 
solution evolves on a smooth manifold has received considerable attention in 
the recent literature [150,161,170]. 
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7 Physical applications 



Prom previous sections it should be clear that ME has a strong bearing on both 
Classical and Quantum Mechanics. As far as Classical Mechanics is concerned 
this has been most explicitly shown in section 3.4. On its turn Quantum 
Mechanics has been repeatedly invoked as a source of applications, among 
others, in sections 2.9 and 4. In this section we present in a very schematic 
way and with no aim at completeness some applications of ME in different 
areas of the physical sciences. This will show that over the years ME has been 
one of the preferred options to deal with equation (4) which, under different 
appearances, pervades the entire field of Physics. In the works mentioned 
here, almost exclusively analytical methods are used and, in general, one must 
recognize that in most, if not all, cases listed only the first two orders of the 
expansion have been considered. In very specially simple applications, due to 
particular algebraic properties of the operators involved, this happens to be 
exact. Only with the more recent advent of the numerical applications, as has 
been emphasized in sections 5 and 6, has the expansion been carried in a more 
systematic way to higher orders. 

7. 1 Nuclear, atomic and molecular physics 

As far as we know the first physical application of ME dates back to 1963. 
Robinson [203] published a brand new formalism to investigate multiple Cou- 
lomb excitations of deformed nuclei. As a matter of fact, he states explicitly 
that only after completion of his work he discovered the ME. His derivation 
of ME formulas is certainly worth reading. 

The Coulomb excitation process yields information about the low lying nuclear 
states. Prior to Robinson work, the theory was essentially based on pertur- 
bation expansions which requires that the bombarding energy is kept so low 
that no nuclear reaction takes place. Even worse, if heavier ions are used as 
projectiles the electric field exerted on the target nucleus is so strong that 
perturbation methods fail. 

The work by Robinson improved the so-called at that time sudden approxi- 
mation^ which is equivalent to the assumption that all nuclear energy levels 
are degenerate. Results are reported in that reference for rotational and vi- 
brational nuclei. 

As representatives of the applications of ME in the field of Atomic Physics we 
mention several types of atomic collisions. The ME is used in [79] to derive 
the transition amplitude and the cross section for K-shell ionization of atoms 
by heavy-ion impact. This is an important process in heavy-ion physics. The 
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theoretical investigations of these reactions always assumed that the projectile 
is a relatively light ion such as a proton or an a particle. The use of ME allowed 
to extend the studies to the ionization of light target atoms by much heavier 
projectile ions. 

In [240,241] the ME is applied to study the time-evolution of rotationally 
induced inner-shell excitation in atomic collisions. In this context the intcrnu- 
clear motion can be treated classically and the remaining quantum-mechanical 
problem for the electronic motion is then time-dependent. In particular, in 
[240] exphcit results for Ne+Ne coUisions are given as well as a study of the 
convergence properties of ME with respect to the impact parameter. 

The ME is applied in [109] to the theoretical study of electron-atom collisions, 
involving many channels coupled by strong, long-range forces. Then, as a test 
case, the theory is applied to electron-impact excitation of the resonance tran- 
sitions of Li, Na and K. Computations up to second order are carried out and 
the cross sections found are in good agreement with experimental data for the 
intermediate-energy range. 

The following examples illustrate the use of ME in Molecular Physics. In 
[38] it is applied for the first time to the theory of the pressure broadening 
of rotational spectra. Unlike the previous approaches to the problem, the S- 
matrix obtained is unitary. As a consequence of it the relative contributions on 
the linewidth of the attractive and repulsive anisotropy terms in the interaction 
potential may be calculated. 

Floquet theory is applied in [169] to systems periodic in time in the semi- 
classical approximation of the radiation-quantum-molecule interaction in an 
intense field. The paper contains an interesting discussion about the appro- 
priateness of the Schrodinger and Interaction pictures. One and two-photon 
probability transitions are obtained up to second order in ME. Noteworthy, 
formulas through fifth order in ME are given, in a less symmetrical form. 

In [210] it is explored the applicability of ME to the multiphoton excitation 
of a sparse level system for which the rotating wave function approximation is 
not applicable. This reference provides a method of treating the time-evolution 
of a pumped molecular system in the low energy region, which is characterized 
by a sparse distribution of bound vibrational states. 

1.2 Nuclear magnetic resonance: Average Hamiltonian theory 

By far this is the field where ME has been most systematically used and so 
we consider it apart. From elementary quantum mechanics it is known that a 
constant magnetic field breaks the degeneracy of the energy levels of an atomic 
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nucleus with spin. If the nuclear spin is s then 2s + 1 sublcvcls appear. In a 
sample these states are occupied according to Boltzman distribution with a 
exponentially distributed population. When a time-dependent radio-frequency 
electromagnetic field of appropriate frequency is applied then energy can be 
absorbed by certain nuclei which are consequently promoted to higher levels. 
This is the physical phenomenon of Nuclear Magnetic Resonance (NMR). 

It was Evans [83] and Haeberlen and Waugh [98] who first applied the ME to 
NMR. Since that time, the ME has been instrumental in the development of 
improved techniques in NMR spectroscopy [36] . 

The major advantage of NMR is the possibility of modifying the nuclear spin 
Hamiltonian almost at will and to adapt it to the needs of the problem to be 
solved [82] . This manipulation requires an external perturbation of the system 
that can be either time-independent (changes of temperature, pressure, sol- 
vents, etc.) or time-dependent (sample spinning, pulsed radio-frequency fields). 
In the later context, the concept of average Hamiltonian provides an elegant 
description of the effects of a time-dependent perturbation applied to the sys- 
tem. It was originally introduced in NMR by Waugh [82,234] to explain the 
effects of multiple-pulse sequences. 

The basic idea of average Hamiltonian theory, for a system governed by H{t), 
consists in describing the effective evolution within a fixed time interval by an 
average Hamiltonian H. The theory states that this is always possible provided 
H{t) is periodic. The average Hamiltonian depends, however, on the beginning 
and the end of the time interval observed. It is right the average Hamiltonian 
H which is obtained by means of the ME. 

When the total Hamiltonian splits in a time-independent and a time-dependent 
piece, H(t) = Hq + Hi{t), with Hi{t) periodic, an interesting new picture 
is used, labeled toggling frame. It certainly reminds the Interaction Picture 
defined in equation (91) but is rather different. In (89) the operator G{t) 
associated to the toggling frame is given by the time-ordered expression 



and the key point here is whether the formal time-ordering is solvable. 

As already mentioned the interplay between NMR and ME has been fruitful 
along the years and acted in both directions. To prove that it is still ahve we 
quote two recent papers directly dealing with that mutual interaction. In [228] 
the relevance of ME through NMR for the new field of quantum information 
processing and computing is envisaged. The authors of [229] have recently 
explored the fourth and sixth order of ME to design a software package for the 
simulation of NMR experiments. Although their results are not yet conclusive 




(355) 
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their work shows the vitahty of the ME. 



1.3 Quantum Field Theory and High Energy Physics 

The starting point of any quantum field theory (QFT) calculation is again 
equation (4) which is conventionally treated by time-dependent perturbation 
theory. So the first question which arises is the connection between ME and 
Dyson-type series. This has already been dealt with in subsection 2.4. The 
main advantage of the first one is, as has already repeatedly pointed out, that 
the unitary character of the evolution operator is preserved at all orders of 
approximation. In the historical development of QFT it was, however, Dyson 
approach what was followed. The lost of unitarity was not thought to be of 
great relevance in front of the problems presented by the infinities appearing 
all over the place. Once renormalization idea was introduced, this awful aspect 
of the theory was put also under control. The results were, from the point of 
view of the calculation of observable magnitudes, an unprecedented success: 
the agreement between experimental results and their theoretical counterparts 
was impressive. 

So no wonder if alternatives to Dyson series, such as ME, did not see popular 
acceptance. However, during the years there has been interesting developments 
involving ME in the context of field theory. In particular its use has shown 
to imply a re-ordering of terms in the calculations in such a way that some 
infinities do not appear and so make not necessary the introduction of coun- 
terterms in the Hamiltonian. This is what happens for example in [217] where 
models are built in which ultraviolet divergence appear neither in the Hamil- 
tonian nor in the S-matrix. In principle the results are valid for relativistic 
field theories with any particle content and with minimal assumptions about 
the form of the interaction. 

ME as an alternative to conventional perturbation theory for quantum fields 
has also been studied in [56] where normal products. Wick theorem and the 
like are used to deduce graphical rules d la Feynman for the terms fti for 
any value of i. This has been proved [58] helpful in the treatment of infrared 
divergences for some QED processes such as the scattering of an electron on 
an external potential or the bremsstrahlung of one hard photon, both cases 
accompanied by the emission of an arbitrary number of soft photons. An inter- 
esting feature of the ME based approach is that the theory is free form infrared 
and mass divergences as a consequence of the unitary character of the approx- 
imate time-evolution operator [58,59]. The method is simpler than previous 
techniques based on re-summation of the perturbation series to get rid of those 
divergences. Furthermore, in contrast with the usual treatment, the resolution 
of the detector is not an infrared regularization parameter. An application to 
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Bhabha scattering (clastic electron-positron scattering) is developed in [57]. 
The difficulties of extending the results to Quantum Chromodynamics are 
commented in [56]. 

Recently, an extension of the Magnus expansion has also been used in the 
context of Connes-Kreimer's Hopf algebra approach to perturbative renor- 
malization of quantum field theory [54,55]. In particular, in [78], it is shown 
that this generalized ME allows one to solve the Bogoliubov-Atkinson recur- 
sion in this setting. 

In the field of high energy physics ME has also found applications. Next we 
just quote two instances: one referring to heavy ion collisions and the other to 
elementary particle physics. 

In collision problems the unitarity of the time evolution operator imposes some 
bound on the experimentally observable cross sections. When these magni- 
tudes are theoretically calculated one usually keeps only the lowest orders in 
conventional perturbation theory. This may be harmless at relatively low en- 
ergies but it may lead to unitarity bounds violation as the energy increases. 
The use of a manifestly unitary approximation scheme is then necessary. ME 
provides such an scheme. In heavy ion collision at sufficiently high energy 
and given kinematic configuration (small impact parameter) that violation 
is produced for e+e~when analyzed in the lowest-order time-dependent per- 
turbation theory. In [111] a remedy for this situation was advanced by the 
use of first order ME. It is discussed how most theoretical approaches are 
based on either lowest-order time- dependent perturbation theory or the Fermi- 
Weizsaker- Williams method of virtual photons. These approaches violate uni- 
tarity bounds for sufficiently high collision energies and thus the probability 
for single-pair creation exceeds unity. With some additional assumptions a 
restricted class of diagrams associated with electron-positron loops can be 
summed to infinite order in the external charge. The electron-positron transi- 
tions amplitudes and production probabilities obtained are manifestly unitary 
and gauge invariant. 

In recent years there has been a great interest in neutrino oscillations and its 
closely related solar neutrino problem. The known three families of neutrinos 
with different flavors (electron, muon and tau) were experimentally shown 
to be able of converting into each other. The experiments were carried out 
with neutrinos of different origins: solar, atmospheric, produced in nuclear 
reactors and in particle accelerators. Here oscillation means that neutrinos 
of a given flavour can, after propagation, change their flavour. The accepted 
explanation for this phenomenon is that neutrinos with a definite flavour have 
not a deflnite mass, and the other way around. Let us denote by \va) the 
neutrinos of deflnite flavour with a the flavor index (i.e., electron, muon, tau) 
and by the neutrinos with well defined distinct masses mj, i — 1,2,3 . 
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Then the previous assertion means that \i'a) will be a linear combination of 
the different . 

As neutrinos with different masses propagate with different velocities this 
mixing allows for flavour conversion, i.e. for neutrinos oscillations. ME enters 
the game in the solution of the evolution operator in one basis. If one neutrino 
"decouples" from the other two then the problem reduces to one with only 
two effective generations. Mathematically it is similar to the two level system 
studied in Section 3. The reader is referred to [67,68,69,222] for details of the 
calculations for two and three generations. 



7.4 Electromagnetism 



The Maxwell equations govern the evolution of electromagnetic waves. The 
equations are linear in its usual form and, although they have been extensively 
studied, the complexity to obtain approximate solutions is rather significant. 
When reformulating the equations for a given problem, where the geometry, 
the boundary conditions, etc. are considered and appropriate discretisations 
are taken into account, it is frequent to end with linear non-autonomous equa- 
tions, so that the Magnus expansion can be of interest here. 

To illustrate some possible apphcations, let us consider the Maxwell equations 

1 _ ^ 
= __VxE 

-^-VxH--JW 

where H, E, J arc the magnetic and electric field intensities, and the current 
density, respectively, /i is the permeability and e is the permittivity. After 
space discritisation, these equations turn into a large linear system of non- 
homogeneous equations and Magnus integrators can in principle be applied. 
Time dependent contributions can also appear from boundary conditions or 
external interactions. In some cases J = o"E [30,126], so that, if the conduc- 
tivity a is not constant, Magnus integrators can be useful. 

Let us now consider the frequency domain Maxwell equations (with J = 0) 

V X E = iwuH 

(357) 

V X H = -iweE 

where w is the angular frequency. These equations are of interest for time- 
harmonic lightwaves propagating in a wave-guiding structure composed of 
linear isotropic materials. If one is interested in the x and y components of 



145 



H and E, and how they propagate in the z direction, the equations to solve, 
after appropriate discretisation, take the form [154] 



(iu dw 

— iwe— = A\z)n ^ —iwa— = B(z)u. (358) 
dz dz 

Here u, v are vectors and A, B matrices depending on z, which in this case play 
the role of evolution parameter. A fourth-order Magnus integrator has been 

used in [154]. From section 5 we observe that higher order Magnus integrators, 
combined with splitting methods could also lead to efficient algorithms to 
obtain accurate numerical results with preserved qualitative properties of the 
exact solution. 



7. 5 Optics 



In the review paper [60] one can find references to some early applications of 
ME to Optics. For example to Hamiltonians involving the generators of SU(2), 
SU(1, 1) and Heisenberg-Weyl groups with applications to laser-plasma scat- 
tering and pulse propagation in free-electron lasers. Here as representatives of 
the more modern interest of ME in Optics we quote two applications referring 
to Helmholtz equation and to the study of Stokes parameters. 

Helmholtz equation in one spatial dimension with a variable refractive index 
n{x) reads 

iP"{x) + k^n'^{x)^{x) = 0, (359) 
where k is the wavenumber in vacuum. 

Recently this time-honored wave equation has been treated in two different 
ways, both using ME. Prom a more formal point of view, in [134] Helmholtz 
equation is analyzed following the well known procedure followed by Feshbach 
and Villars to convert the second order relativistic quantum Klein-Gordon 
differential equation for spin-0 particles in a first order differential equation 
involving two components wave functions (the original wave function and its 
time derivative) . The evolution operator for Helmholtz equation is then a 2 x 2 
matrix which evolves according to the fundamental equation (4) with the only 
difference that now the evolution parameter is x instead of t. In [134] the whole 
procedure is explained and the main physical consequence, which amounts to 
the addition of correcting terms to the Hamiltonian, is discussed in the case 
of an axially symmetric graded-index medium, i.e. one in which the refractive 
index is a polynomial. 

Helmholtz equation has also been investigated with the help of ME in [153,154,155]. 
Here the propagation in a slowly varying waveguide is considered and the 
boundary value problem is converted into an initial value problem by the in- 
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troduction of appropriate operators which are shown to satisfy equation (4). 
Numerical methods to fourth order borrowed from [120] are then used. 

Since mid 19th century the polarization state of light, and in general electro- 
magnetic radiation or any other transverse waves, is described by the so called 
Stokes parameters which constitute a four-dimensional vector S{uj) depending 
on the frequency u. When the light traverses an optical element which acts 
on its polarization state the in and out Stokes vectors are related by 

= M(a;)Sin(u;), (360) 

where the 4x4 matrix M{uj) is called the Mueller matrix. It can be proved 
[201,202] that it satisfies the equation 

M'{uj) = H{uj)M{u), M{uJo) = Mo, (361) 

where now the prime denotes derivative with respect to the real independent 
variable uo. For systems with zero polarization-dependent loss (PDL) and no 
polarization mode dispersion (PMD) H{uj) is constant whereas with PDL and 
PDM the previous equation is just our equation (4) and the appropriateness 
of ME is apparent. The matrix H(u}) in this application has an Hermitian 
and a non-Hermitian component. ME has allowed a recursive calculation of 
successive orders of the frequency variation of the Mueller matrix. This yields 
PMD and PDL compensators that counteract the effects of PMD and PDL 
with increased accuracy. 

Also related to the use of Stokes vector one can mention the so-called radiative 
transfer equation for polarized light. It is relevant in Astrophysics to measure 
the magnetic fields in the Sun and stars. That equation gives the variation of 
S{z) with the light path z 

±S{z) = -K{z)S{z)+J, 

where S is the Stokes vector, 7^ is a 4 x 4 matrix which describes absorp- 
tion in the presence of Zeeman effect and J stands for the emission term. In 
[148,149,213] ME is used to obtain an exponential solution. 

7.6 General Relativity 

To illustrate once more the pervasive presence of the linear differential equa- 
tion (4) let us mention reference [168] in which the aim is to determine the 
time elapsed between two events when the space-time is treated as in General 
Relativity. Then it turns out to be necessary to solve a two-point boundary 
value problem for null geodesies. In so doing one needs to know a Jacobian 
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whose expression involves a 8 x 8 matrix function obeying the basic equation 
(4). In [168] an eighth order numerical method from [22] is used, which is 
proved to be an efficient scheme. 



7.7 Search of Periodic Orbits 

The search of periodic orbits for some non-linear differential autonomous equa- 
tions, x' = f(x), X e R*^ is of interest in Celestial Mechanics (periodic orbits 
of the N-body problem) as well as in the general theory of dynamical sys- 
tems. Due to the complexity of this process, it is important to have efficient 
numerical algorithms. 

The Lindstedt-Poincare technique is frequently used to calculate periodic or- 
bits. An iterative process proposed in [231] consists in starting with a guessed 
periodic orbit, and this guess is subsequently improved by solving a correla- 
tion non-autonomous linear differential equation. The numerical integration 
of this equation is carried out by means of Magnus integrators. 



7.8 Geometric control of mechanical systems 

Many mechanical systems studied in control theory can be modeled by an 
ordinary differential equation of the form [33,132] 

m 

x'(i) = fo(x(t)) + Y.<t) f,(x(i)), (362) 

i=l 

initiahzed at x(0) = p. Here x e R'' represents all the possible states of the 
system, fj are (real) analytic vector fields and the function u = (-Ui, . . . ^Um) 
(the controls) are assumed to be integrablc with respect to time and taking 
values in a compact subset U C M™. The vector field fo is called the drift 
vector field, whereas fj, i > 1, are referred to as the control vector fields. 
When fo = 0, the system (362) is called 'without drift', and its analysis is 
typically easier. 

For a given set of controls {ui}, equation (362) with initial value x(0) is nothing 
but a dynamical system, which can be analyzed and (approximately) solved 
by standard techniques. In control theory, however, one is interested typically 
in the inverse problem: given a target x(T), find controls {ui} that steer from 
x(0) to x(T) [132], perhaps by following a prescribed path. Just to illustrate 
these abstract considerations, a typical problem could be the determination 
of a set of controls that drive the actions of a robot during a task. 
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The first step is to guarantee that there exists a solution. This is the problem 
of controllability. To characterize the controllability of linear systems of the 
form 

x.'{t) = A(t)x(t) + B{t)u{t) (363) 

is a relatively simple task thanks to an algebraic criterion known as the Kalman 
rank condition [33,129]. This issue, however, is much more involved for the 
nonlinear system (362) [132]. 

The interest of ME in control theory, as it has been already discussed in 
section 3.5, stems from the approximate ansatz it provides connecting the 
states x(0) and x(T). Thus, the Magnus expansion can be used either to 
predict a state x(T) for a given control u or to find reasonable controls u 
which made reachable the target x(T) from x(0). Of course, many sets of 
controls may exist and it raises questions concerning the cost of every scheme 
and consequently the search for the optimal choice. For instance, the ME has 
been used in non-holonomic motion planning of systems without drift [72,73]. 
Among non-holonomic systems there are free-floating robots, mobile robots 
and underwater vehicles [73,185]. 

In the particular case of linear quadratic optimal control problems (appear- 
ing in in engineering problems as well as in differential games) a given cost 
functional has to achieve a minimum. When this happens, eq. (363) can be 
written as [80,200] 

x' = M{t,K{t))x, (364) 

where K{t) e R'^^'^ has to solve a Riccati differential equation similar to 
(339) with final condition K(T) — Kf. In other words, the Riccati equation 
has to be integrated backward in time and then to use it as an input in 
(364). As mentioned, the Riccati differential matrix equation has received 
much attention [24,63,80,127,133,200,211], but an efficient implementation to 
this problem requires further investigation, and methods from ME can play 
an important role. 



8 Conclusions 

In this report we have thoroughly reviewed the abiding work on Magnus expan- 
sion carried out during more than fifty years from very different perspectives. 

As a result of a real interdisciplinary activity some aspects of the original 
formulation have been refined. This applies for example to the convergence 
properties of the expansion which have been much sharpened. 

In other features much practical progress has been made. This is the case of 
the calculation of the terms of the series both explicitly or recurrently. New 
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techniques, like the ones borrowed from graph theory, have also profitably 
entered the play. 

Although originally formulated for linear systems of ordinary differential equa- 
tions, the domain of usage of ME has enlarged to include other types of prob- 
lems with differential equations: stochastic equations, nonlinear equations or 
Sturm- Liouville problems. 

In parallel with this developments in the mathematical structure of the ME 
the realm of its applications has also widen with the years. It is worth stressing 
in this respect the versatility of the expansion to cope with new applications 
in old fields, like NMR for instance, and at the same time its capability to gen- 
erate new contributions, like the generation of efficient numerical algorithms 
for geometric integrators. 

All these facts, historical and present, presented and discussed in this report 
strongly support the idea that ME can be a very useful tool for physicists. 
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